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STANISŁAW SKOCZOWSKI 


IDENTYFIKACJA I MODELOWANIE UPROSZCZONE 
OBIEKTÓW ELEKTROTERMICZNYCH NA PODSTAWIE 
POCZĄTKOWEJ FAZY CHARAKTERYSTYKI 
ROZGRZEWU 


Instytut Automatyki Przemysłowej Politechniki Szczecińskiej, 
ul. Gen. Sikorskiego 37, 70-313 Szczecin 


Streszczenie: Omówiono podstawy teoretyczne oraz wyniki badań eks- 
perymentalnych nowej metody deterministycznej identyfikacji modeli uproszczonych 
w zastosowaniu do pieców rezystancyjnych. 


Słowa kluczowe: modelowanie — identyfikacja — elektrotermia. 


1. WPROWADZENIE 


Celem regulacji temperatury w piecach oporowych posługujemy się najprost- 
szymi liniowymi modelami o stałych skupionych, najczęściej w postaci inercji 
pierwszego rzędu z opóźnieniem. To Świadomie czynione uproszczenie rzeczywi- 
stych procesów elektrotermicznych, które pozbawione działania różniczkujące- 
go są z natury nieliniowymi o zależnych od temperatury parametrach [4—6] 
układami o stałych rozłożonych, może okazać się przyczyną słabej jakości 
regulacji temperatury, niewystarczającej w wielu procesach technologicznych. 

W elektrycznych piecach oporowych występująca na ogół nieliniowość 
statyczna o charakterze degresywnym ujawnia się [13, 14] różnymi stałymi 
czasowymi dla załączenia i wyłączenia mocy pieca. Zjawisko różnych stałych 
czasowych dla nagrzewania i stygnięcia, znane w elektrotermii od dawna [2, 7, 
9—11], zostało także zauważone w innych dziedzinach [3, 8]. 

W pracach [10, 11] wykazano, że do analizy regulacji on-ofj z uwzględ- 
nieniem przeregulowania, musi być użyty model obiektu elektrotermicznego 
co najmniej drugiego rzędu. Zalety modelu n-tego rzędu o jednakowej stałej 
czasowej, analizując w kategoriach procesów stochastycznych redukcję wa- 
riancji na wyjściu, przedstawiono w pracy [12]. Ostatnio [1] stwierdza się, 
że model pierwszego rzędu z opóźnieniem nie jest reprezentatywny dla 
większości procesów przemysłowych, które powinny być modelowane wyż- 
szym rzędem. 


Rozpatrując proces elektrotermiczny nagrzewania wsadu z regulacją tem- 
peratury, dostrzega się wyraźnie więcej niż jeden magazyn energii, możemy 
bowiem wyróżnić co najmniej trzy podukłady: czujnik temperatury. w obudo- 
wie, wsad i komora pieca, z których każdy mógłby być modelowany inercją 
z opóźnieniem. Tak więc pojawia się problem identyfikacji rzędu przyjętego 
prostego modelu transmitancyjnego wieloinercyjnego. Prezentowana tu meto- 
da (15—20) opiera się na próbkowaniu dyskretnym początkowej fazy od- 
powiedzi skokowej rzeczywistego procesu przy założeniu znajomości czystego 
opóźnienia. W przypadku pieców oporowych może ona okazać się szczególnie 
przydatna, ponieważ zwykle dla każdego pieca rejestruje się rutynowo, zgodnie 
z normą, tzw. charakterystykę rozgrzewu, która w rzeczywistości jest obszer- 
nym fragmentem odpowiedzi skokowej obiektu. Umożliwiło to autorowi 
weryfikację praktyczną badań teoretycznych. Ponadto metoda oparta na 
pomiarze próbek odpowiedzi położonych blisko siebie jest odporna na wpływ 
nieliniowości statycznej pieca. Zaletą metody jest również szybkość uzyskania 
wyników i jej prostota. Największy wpływ na dokładność wyników mają 
zakłócenia losowe nałożone na próbkowany sygnał użyteczny. Wykorzystanie 
pojęcia Średniej stałej czasowej [23] przy znajomości rzędu pozwala, dla 
procesów aperiodycznych stabilnych, na przejście od modelu rzędu n o jed- 
nakowej średniej stałej czasowej do modelu wieloinercyjnego rzędu n o różnych 
stałych czasowych z określonym ich rozrzutem [18]. Dla przybliżonego 
modelu liniowego wyższego rzędu można zaprojektować regulator PID 
odporny na zmiany parametrów modelowanego procesu, projektując nastawy 
regulatora ze względu na margines fazy i amplitudy [21]. 

Próby praktyczne wykorzystania prezentowanej metody do modelowania 
rzeczywistych procesów elektrotermicznych dały wyniki pozytywne [16—18, 
20, 22]. 


2. PODSTAWY TEORETYCZNE KONSTRUKCJI 
UPROSZCZONYCH MODELI 


2.1. MODELE TRANSMITANCYJNE PROCESÓW APERIODYCZNYCH 


Przy założeniu liniowości proces można scharakteryzować transmitancją: 


b„S”+ ... +b, s+b 
G(s) = - m —— 1 0 
$Fd, yk «+ Pdy S$ FA 


eT s", (1) 


Definiuje się pojęcie Średniej stałej czasowej T, która spełnia relację: 


l _ Pit: +P) _ —(6G17-:- +0, 


T n n (2) 


gdzie o, = Retp;y, p, są biegunami transmitancji (1). 


Przy tym ważne są zależności: 


pr n 
T=—=——, (3) 
dh-1 l 
is FB 
Tmin < T < dyać (4) 


Definiuje się też współczynnik rozrzutu A w postaci: 


A=n-1]7=%, 0</<l, (5) 


max 


co oznacza, że poszczególne stałe czasowe są „rozpięte” pomiędzy Ty» I Tmax 
wg zależności: 


= WEŃ, os Bt. (6) 


Zatem stabilny proces aperiodyczny (m = 0) może być zastąpiony modelem 


S(s) = (7) 


(1+sT)Y" 


a przy znajomości rozrzutu A można przejść do dokładniej aproksymującego 
proces modelu 


M(s)==-———. (8) 
II G+sT) 


Zauważmy, że w przypadku znanego opóźnienia r modele (7) i (8) mogą być 
prosto skorygowane przez uwzględnienie dodatkowego członu e ”. 


2.2. SZACOWANIE PARAMETRÓW MODELI 
Odpowiedź skokowa modelu (1) jest odwrotną transformatą Laplace'a: 
G (s 
hQ="L" zł (9) 


Dzieląc licznik wyrażenia (1) przez mianownik i dokonując operacji (9), 
przy założeniu c =0, uzyskuje się szereg 


tę r+1l 
Mt) =Pb => 
(0 = bajt bn-1—bn a, Vy 
pro 
Ebm 2 bn 1 01-1 Fb (Gs 1— 9-2) "jp 20 (10) 


gdzie r = n—m. 


10 


Szereg ten jest sumą parabol i dla niewielkich wartości argumentu t jest 
silnie zbieżny. Biorąc stosunek dwu wartości dyskretnych odpowiedzi skoko- 
wej dla t, i t,, uzyskuje się dla niewielkich relatywnie wielkości t;/T, t, < t, 


zależność: 
ka (el K. (11) 
ADOBE 


gdzie funkcja $© jest w przybliżeniu równa jedności dla niewielkich relatywnie 
okresów próbkowania. Powinno być spełnione t, < 0,7T. Logarytmując zwią- 
zek (11) i ograniczając się do modeli (7) i (8), to jest dla m =0, r = n, można 
napisać 


ZE ą | (12) 


gdzie ó, nazwiemy „błędem rzędu”, a n, jest określone równaniem: 
h (t, 
In h A 
z ŻE, (13) 
J 


hy = 
t 
In = 
L, 

Rząd n jest zaokrągleniem w górę n, do liczby całkowitej: 
n =Ent nv). (14) 
Znając rząd n, można w przybliżeniu na podstawie zależności (10) i (11) 

oszacować pozostałe parametry modelu (7): 


n 
n+1 t,| h(t.) | A 
„(GVE60_, Le kle) 
(eJra"" 
-_h*(t) RUE 
Fm (a) (6 . 


„Srednie wzmocnienie” określone związkiem (16) może posłużyć do oszaco- 
wania wzmocnienia k. Ogólnie biorąc, w transmitancji (1) wzmocnienie 


Tx (15) 


k = b,/aq. Dla modelu (7) a, = 1/] | T;, co uwzględniając i biorąc pod uwagę 
i=1 
wyrażenia (3) i (6), otrzymujemy 


R 
ROCZEK a 


T" TYPEZEWEJCEJOWE (17) 
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Ze względu na zakłócenia, które mogą być nałożone na próbkowaną 
odpowiedź skokową, można wykorzystać całą serię N próbek, z której 
uzyskuje się N(N —1)/2 par t;, t,, a zatem i tyleż wyników oszacowań, które po 
odrzuceniu śalnsch mogą zostać uśrednione (odnosi się to do ny, Ti k). 
Także sygnały przed próbkowaniem mogą być filtrowane. 

Funkcja $ może być z dużą dokładnością aproksymowana [19] dla modeli 
(7) i (8) zależnością: 

nh żę 


bze PT, (18) 


co uwzględniając w równaniu (12), daje 


O =n=-n, Z ki ć (19) 


ponieważ 0, może być łatwo określone podczas identyfikacji rzędu, to wzór ten 
może posłużyć także do oszacowania T, niezależnie od wzoru (15). Do 
oszacowania parametrów modelu (8) niezbędne jest oszacowanie współczyn- 
nika rozrzutu A oraz dodatkowo niezbędna jest znajomość jednej ze stałych 
czasu T;, np. T;„ lub T,,,. Często w praktyce wymóg ten może być spełniony, 
bo przykładowo może być znana stała czasu czujnika temperatury zain- 
stalowanego w komorze obiektu elektrotermicznego stanowiącego razem 
z czujnikiem obiekt identyfikacji. 

Opierając się na zależnościach (3) i (5) można napisać równania dla 
obliczenia A przy znajomości albo T,,„, albo Ta, [18, 19]: 


T, 
n1 4, 1+1—n=F=0, 
A „+ŁA+ n T % 
4 1 paz |, +4+1=0 
T Ka =0. 


Na podstawie układu (20), dla 4 >0, T,,,, > Tjn otrzymujemy oszacowanie: 


= PAS = T. (21) 


s|> 


Dysponując dodatkową informacją o czasie t,, w którym odpowiedź 
skokowa w punkcie przegięcia osiąga maksymalną pochodną (dh/dt),„,, można 
wykonując obliczenia numeryczne wyznaczyć 4. Przykładowo, dla n = 2 będzie 
to równanie 


„żuk 
R 1ESAFA 


(22) 


a 


Podobne wyrażenia dla n=3 i n=4 podano w pracy [19]. 


2.3. MODEL DYSKRETNY 


Uwzględniając istnienie ekstrapolatora zerowego rzędu w cyfrowym sys- 
temie sterowania oraz okres próbkowania 4, model dyskretny może być 
przedstawiony z zastosowaniem przekształcenia Z w postaci impulsowej 
funkcji przejścia 


H()=Z je o Ć | = Gspz * (23) 
gdzie 
1 
d=2, (24) 
A(z')=l+a,z "+... +0,z ”, (25) 
B(z79)=P,z71+...+$,z"", (26) 
z ! =exp(—45). (27) 


Dobór okresu próbkowania 4 jest problemem oddzielnym. Powinien on 
być dobrany tak, by spełniona była nierówność 


4 < 0,3 dniaj 


ale nie może być zbyt mały, gdyż prowadzi to do degeneracji modelu (23). 
Współczynniki wielomianu A(z ") mogą być w przybliżeniu określone [23] 
na podstawie zależności: 


IA 

a Zryw (-F) KICI, (28) 
Tak więc dla wyznaczenia współczynników wielomianu A(z *) należy znać 
jedynie rząd n i Średnią stałą czasową T. 

Oszacowanie współczynników wielomianu B(z” *) jest bardziej złożone. Na 
podstawie wyrażenia (28) i przytoczonych w pracy [23] zależności obliczono tu 
dla modelu (7) S(s) współczynniki wielomianów A(z ') i B(z !) modelu 
dyskretnego, dla rzędu n= 1, 2, 3 i 4. Przedstawiono je w tabeli 1. 

Należy podkreślić, że dla obliczenia parametrów modelu dyskretnego 
niezbędna jest identyfikacja n, Ti by. Współczynnik by w transmitancji (1) 
może być określony na podstawie zależności (16): 


Kx h* (t.) Gz. (29) 


h(t) At) ti 


Dysponując modelem dyskretnym, możliwe jest także inne szacowanie 
wzmocnienia k. Współczynnik a, w mianowniku transmitancji (1) może być 
oszacowany [23] z użyciem modelu (7): 

A(1) 


a, AE(2) (30) 
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Tabela 1. Współczynniki modelu dyskretnego dla rzędu 1, 2, 314 
Współczynniki A(z ") Współczynniki BZ”) 


1 


Ks 
Il 
GG 

[>] 

NI 
TEn 
— 
| 
s RZ 
NIE 


2T 
LA. 24 
| Poz OaT 
SĄ! AŻ 1-24 > 
Ha oka 
A3/ 34 
p =bof :7) 
3 
114? : 
; a sr) 
1214*/, 44 8 4 
armory (-57) 
_, 1214*/,_44 
(3 4 WaS 
gdzie 
A()=l+a, +... +a,, (31) 
a 
5 n4 
E(AJZuzr=czp| 2]. 3 
(4) 0 ap 7)) (32) 
Wtedy wzmocnienie 
- b 
kz—. (33) 


Oczywiście znajomość rzutu 4 prowadzi do dokładniejszego oszacowania 
ay, a co za tym idzie, do lepszego oszacowania wzmocnienia — wzór (17). 
Należy podkreślić, że identyfikowany model dyskretny jest budowany w celu 
sterowania i regulacji temperatury, stąd przyjęty okres próbkowania 4 będzie 
zwykle większy od okresu próbkowania w czasie identyfikacji. 


2.4. PRZYBLIŻONE SZACOWANIE NIELINIOWOŚCI STATYCZNEJ 


Rozważając zagadnienie modelowania nieliniowości statycznej obiektu 
elektrotermicznego opieramy się na prostym schemacie przedstawionym na 
rycinie 1. W czasie identyfikacji obiektu rzeczywistego dostępne są pomiarowo 
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Liniowy Model X(1) 
model charakterystyki p—— 
dynamiki nieliniowej 


Ryc. 1. Model obiektu z uwzględnieniem nieliniowości 


sygnały U (t) i X (t). Zauważmy ponadto, że w przypadku pieców oporowych 
zwykle nie można eksperymentalnie przez pomiary bezpośrednie wyznaczyć 
temperatury ustalonej pieca, odpowiadającej załączeniu na stałe mocy znamio- 
nowej, gdyż groziłoby to zniszczeniem pieca. Specyfiką tych obiektów jest także 
to, że można zwykle załączyć skokowo — w sposób idealny — pełną moc 
znamionową lub jej część na wystarczająco długi okres. W pracach [11, 13, 14] 
wykazano, że stosując model (7) przy 1>0, n=l|1 (model Kiipfmiilera), 
w dużym przybliżeniu za model charakterystyki statycznej obiektów elektro- 
termicznych można przyjąć równanie: 
u 

— Ati(l=n) e) 
gdzie x i u są względnymi wartościami średnimi ustalonymi, a tzw. stopień 
nieliniowości „u jest określany wzorem: 


(35) 


w którym «, i «„_ są wartościami maksymalnej pochodnej dx (t)/dt = dh(t)/dt 
odpowiedzi skokowej odpowiednio dla stanu nagrzewania i stygnięcia. Wspo- 
mniano o tym wcześniej, definiując czas t, osiągnięcia punktu przegięcia 
odpowiedzi skokowej. 

Równanie (34), obowiązujące w stanie ustalonym, dotyczy względnych 
wartości średnich na wyjściu i wejściu, odniesionych do wartości maksymal- 
nych. Podczas regulacji temperatury piec będzie sterowany w stanie ustalonym 
średnią mocą P. Odnosząc P do mocy znamionowej P,,, można napisać 


Ę 


0<= 
Pzn 


1 (36) 


a ponieważ ii = P/P,,, więc nierówność (36) określa zakres zmian u. Sterowa- 
niu w takim zakresie będzie w stanie ustalonym odpowiadać na wyjściu średnia 
temperatura ustalona k. Odnosząc średnią temperaturę ustaloną pieca do 
wzmocnienia k, równego ustalonej temperaturze przy załączeniu na stałe 
(teoretycznie możliwe) pełnej mocy znamionowej, uzyskujemy nierówność 


0 Ea]! (37) 
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Podobnie jak poprzednio, ponieważ x =k/k w stanie ustalonym, więc 
nierówność (37) określa zakres zmian x. Zauważmy, że rozważania w stanie 
ustalonym oznaczają możliwość pominięcia na rycinie 1 dynamiki i uwzględ- 
nienia tylko wzmocnienia k. Zatem przy założeniu, że znany jest stopień 
nieliniowości u określany wzorem (35), można na podstawie równania (34) 
oszacować w przybliżeniu charakterystykę statyczną pieca, a mając choć jeden 
punkt pośredni z pomiarów dla stanu ustalonego (np. wzmocnienie k od- 
powiadające załączeniu połowy mocy znamionowej) można sprawdzić słusz- 
ność tego oszacowania. 

W pracach [13, 14] wykazano, że stopień nieliniowości u spełnia równanie: 


0x 
a 0u| - 
=> RZY = — 2 ą 
Ol, 0x| 8) 
0u|z, 


gdzie x_ odpowiada rzędnej punktu przegięcia odpowiedzi skokowej dla 
stanu stygnięcia, a x, dla stanu nagrzewania. Różniczkując wyrażenie (34) 
otrzymamy 


X Mo 
0u [u+a(l—n)]? 


Zależności (38) 1 (39) mogą posłużyć do kontroli oszacowania charakterystyki 
statycznej. Zauważmy, że podstawowy problem polega na możliwości okreś- 
lenia «„„_ lub x_ na drodze pomiarów, ponieważ praktycznie nie jest możliwe 
zdjęcie pełnej — od k poczynając — charakterystyki skokowej stygnięcia dla 
większości obiektów elektrotermicznych. 

Zauważmy, że w przypadku nieznajomości «„_ dla oszacowania stopnia 
nieliniowości u można [11] posłużyć się przekształconą zależnością (34): 


(39) 


Age. (40) 


Mając choć jeden punkt pośredni pomiarów dla stanu ustalonego, możliwe 
będzie oszacowanie u z zastosowaniem wzoru (40). 


3. PRZYKŁADY ZASTOSOWANIA 


Przedstawione zostaną modele kilku wybranych, rzeczywistych obiektów 
elektrotermicznych z zastosowaniem omówionego sposobu identyfikacji. 

1. Piec elektryczny wgłębny PEGat-950/3 o mocy 90 kW produkcji 
Zakładów Elterma w Świebodzinie. Badanie tzw. czasu rozgrzewu przepro- 
wadzono u wytwórcy zgodnie z PN-73/E-06209. Termopara pomiarowa była 
umieszczona w środku geometrycznym komory pieca. Temperatura w *C była 
odczytana z miernika temperatury regulatora, wyposażonego w swój czujnik 
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temperatury. W tabeli 2 przedstawiono zarejestrowane wyniki z uwzględ- 
nieniem temperatury otoczenia (18 C). Za podstawę obliczeń wzięto wyniki 
pomiarów uzyskanych za pomocą termopary pomiarowej w mV. Wyniki 
obliczeń zestawiono w tabeli 3. Dane pomiarowe podane w tabeli 2 pozwalają 
aproksymować graficznie [19] wartości wzmocnienia (temperatura ustalona 
pomniejszona o wartość temperatury otoczenia) oraz oszacować Ta i t,. 
Ponieważ znana jest temperatura końcowa 3, =932'C dla t, = 240 min, 
można także sprawdzić, w jakim stopniu odpowiedź modelu odbiega od tej 
wartości. Na tej podstawie można przyjąć modele 


1200 


S() % + 518,3)? 4:=-0;1., 
M(s) = 1200 
7 (1+s10)(1 +s100) 


Tabela 2. Dane pomiarowe pieca PEGat-950/3 


Model M (s) daje dla czasu t, = 240 min Jęmoa = 930 C. W tym przypadku, 
przyjmując by = 0,047, uzyskuje się, po przeliczeniu mV na *C, k =405'C, 
a zgodnie z wzorem (17), po uwzględnieniu 4 = 0,1, k = 1227”C. Wzmocnienie 
12007?C przyjęto opierając się na aproksymacji graficznej krzywej rozgrzewu. 
Wartość Tmi;n = 10 spełnia nierówność (21) oraz równanie (20). T,, = 
= (1+4)T/2 dla t, x 25,5 min, A = 0,1 spełnia też równanie (22), a oszacowane 
Imax na podstawie końca zarejestrowanej krzywej rozgrzewu [19] wynosi 
około 90 min. Wyniki te należy uznać za satysfakcjonujące. Model dyskretny 
może być określony po przyjęciu okresu dyskretyzacji 4. Przyjmując 
A=2min, T=18,3 min, b, = 0,053, po przeliczeniu by na ?C, zgodnie 


7 


z wzorami podanymi w tabeli 1 dla n = 2, uzyskuje się model dyskretny pieca 
PEGat 950/3 


2,5442 '+2,359z * 


-1 La J 
H(z )=1—17932-1+0,8042 7 


Z zależności (31) wynika A(1)=0,011, a na podstawie wzoru (32) 
E(4) x 0,89. Zgodnie z oszacowaniem (30) a, £ 3,1:107* i na podstawie 
związku (33) k m 17,2 mV, co daje k m 437?C. Uwzględniając A = 0,1, wzmoc- 
nienie k oszacowane w ten sposób wyniesie x 1322”C. Gdyby do obliczenia 
parametrów modelu dyskretnego wzięto bę = 0,047, niewiele zmieniłyby się 
parametry B w liczniku transmitancji dyskretnej, a tak oszacowane wzmoc- 
nienie KU T2"Ć: 

2. Suszarka elektryczna SDN-IE produkcji Zakładów Elterma z ter- 
moelementem bez obudowy NiCr-Ni umieszczonym w środku geometrycznym 
komory o wymiarach 330 x 300 x 410 mm. 8,, = 250%C, Yptocz = 17,857C. Za- 
rejestrowany w trakcie badania czasu rozgrzewu przebieg daje wartości próbek 
po odjęciu Jptocz JaK W tabeli 4, w której również podano wyniki obliczeń rzędu 
n, wg zależności (13). Zwraca uwagę duża regularność oszacowania ng. Tylko 
7 na 91 wyników powinno być odrzuconych, gdyż są większe od n=2, 
hę = 1,7881. Pełna zdjęta charakterystyka rozgrzewu od Baocz do 9,, pozwala 
na określenie t, = 3900 s, 9, = 244'C, t„ m [1005-1140] s, «,, — 0,083”Cy/s. 
Na podstawie przytoczonych w tabeli 4 t, i h,, stosując zależności (15) 1 (16), 
można przyjąć T = 960 s oraz k = 227?C. Ponieważ t„/T = [1,047-- 1,187], to 
dla t„/T = 1100/960 = 1,146 uzyskuje się A x 0,26, co daje k = 346C. Zatem 
uzyskuje się model przybliżony 

346 


ai 2211 4852026: 
56)  1+:960)7* 


Tabela 4. Dane pomiarowe oraz wyniki obliczeń ny suszarki SDN-1E 


JAOHUDEDUORDLEDE 
|| 0,23 [1,407 [1,518 |1,542 | 1,606 | 1,641 | 1,680 | 1,697 | 1,696 | 1,703 | 1,701 

| Z % 0,61 1,75 |1,788|1,83 1,83 [1,83 |1,82 
Ż 90 | 1,22 1,83 [1,89 1.898 1,873|1,87 |1,855 
4 |100| 1,95 1,98 |2,02 [2,00 |[1,959|1,946|1,919 
5 |150| 3,05 1,947|2,03 |2,01 |1,942|1,927|1,895 
6 | 180 2,139 [2,052 [1,939 | 1,920 [1,823 
7 | 210 1,950|1,816|1,825]|1,790 
8 | 240 1,664 |1,75 |1,723 
9 | 270 1,846 | 1,758 
10 | 300 1,66 
11 | 330 

12 |360 

13 | 390 

14 | 420 


2 — Szczecińskie Roczniki Naukowe t. XI, z. 1 
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3. Suszarka elektryczna SDN-IE produkcji Zakładów Elterma z czuj- 
nikiem temperatury w osłonie. Jest to ta sama suszarka, co omawiana 
w przykładzie 2, a różnica polega na tym, że jednocześnie z rejestrowaniem 
charakterystyki rozgrzewu rejestrowano wskazania cyfrowego miernika tem- 
peratury w regulatorze, wyposażonym w czujnik temperatury Pt100 w osłonie 
ceramicznej i metalowej. W związku z tym temperatury odczytane były 
z niewielką dokładnością — bez miejsca po przecinku. Podano je w tabeli 5, 


Tabela 5. Dane pomiarowe oraz wyniki obliczeń nę MYSIA DEN SDN-1E z czujnikiem w osłonie 


JADBZNZZDEA 


ż 180 2,40 | 2,71 | 2,45 | 2,48 | 2,46 | 2,52 
3 240 j 3,44 | 2,51 | 2,55 | 2,49 | 2,58 
4 270 6 1,46 | 2,02 | 2,11 | 2,30 
5 300 7 2,65 | 2,48 | 2,64 
6 330 9 2,31 | 2,64 
7 360 | 11 3,01 
8 390 | 14 

a 420 | 16 


w której zamieszczono również wyniki obliczeń rzędu n,. Po odrzuceniu 
wyników ny >3i ny <2 ny =2,462+0,1795. Zatem uzyskuje się ewidentny 
wzrost rzędu do n=3, co jest spowodowane zapewne tym, że w trakcie 
identyfikacji został zauważony wpływ inercji czujnika temperatury w osłonie. 
Pomimo ograniczenia do liczb całkowitych wartości próbek h,, uzyskuje się 
dużą powtarzalność wyników ny. Z charakterystyki rozgrzewu rejestrowanej 
także po wyłączeniu mocy przez regulator nastawiony na wartość temperatury 
znamionowej można uzyskać «„_ = 0,0336'C/s. Początek charakterystyki 
rozgrzewu pozwala ocenić «, = 0,09”C/s. Zatem można oszacować stopień 
nieliniowości na podstawie wzoru (35): u = 0,3741. Aby oszacować wzmoc- 
nienie inaczej niż w przykładzie 2, można zastosować metodę, którą uzasadnia 
rycina 2, ukazująca końcowy fragment charakterystyki rozgrzewu. 
Łatwo wykazać, że wzmocnienie k może być oszacowane za pomocą 
zależności: 
ka arów. (41) 


04 0B 
a maksymalna stała czasowa, gdy tylko B jest dostatecznie blisko k, 


k— R 


OB 


(42) 


ESL Ge 


W tym przypadku potrzebne wartości wynoszą: B=244C, ap Z 
= 2,2 C/min, A = 302'C, a, z 3*C/min. Na podstawie związku (41) uzyskuje 
się k =359C, co jest bliskie oszacowaniu w przykładzie 2. Przyjmując tę 
wartość wzmocnienia oraz uzyskane na podstawsie zależności (15) i (19) 
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a 
pk cze 
— 


t tą tj 


Ryc. 2. Ilustracja sposobu szacowania wzmocnienia 


wartości średnie stałej czasowej średniej, odpowiednio T = 670 si T = 600's, 


dające I = 635 s, można zaproponować następujący model suszarki z czuj- 
nikiem w osłonie: 


239 


S(8) X 15635)" 


Zakładając wartość stałej czasowej T,,, = 320 s, na podstawie wzoru (20) 
uzyskuje się A = (0,37, co pozwala na oszacowanie M (5): 


359 


i,  NZO0ZT. 
(1+s320)(1 + s865)(1 + s2340) 


M (5) z 
Aby zbadać, czy powyższy model M (s) daje odpowiednią wartości Dy mod (£,)> 
należy zastosować wzór określający odpowiedź skokową dla n =3. Jest to 
zależność 


t t 
Ram fe TG 
M0 K(M p EOTT"TOTESBZIK* 
PEEBEOB=LE MNSESRSTE 
I 
Bee 
| aa E at (4) 
FWIRSK TE) ( ) 


Dla t, = 3960 s uzyskuje się z równania (43) v, = 237”C. Można przyjąć, że 
TQax suszarki z czujnikiem w osłonie jest bliskie T,,,, suszarki z termoparą 
bez osłony. Zatem po skorygowaniu model zaproponowany w przykładzie 2 
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ma postać: 
359 


——,  N=0,26, 
(1 +s960)* 


S(s) z 


a po uwzględnieniu T,,, = 2330 s, uzyskuje się dla suszarki z termoparą bez 
osłony 


WE 359 
© 7 (14605) (1 + s2330) 
Stosując dla n = 2 wzór 
śę Ra ET o Sęć 
)=kl|1 —ę 1h ź azj, 44 
PA (= Uys | SW 


uzyskuje się temperaturę końcową dla tego modelu bvpmod = 2707C, co można 
uznać za zadowalający wynik, z uwagi na cel identyfikacji. 

4. Piec elektryczny PM-6/1100, Przemysłowy Instytut Elektroniki, 
Warszawa. Zarejestrowane cyfrowo charakterystyki rozgrzewu pieca pustego 
z termoparą w osłonie dla 50% mocy znamionowej i dla 100% mocy 
znamionowej pozwalają zilustrować także metodę szacowania charakterystyki 
statycznej pieca. Okres próbkowania przy rejestracji charakterystyk wynosił 
10 s. Zarejestrowano także charakterystykę stygnięcia pieca. Dla 50% mocy 
znamionowej uzyskano dane (po odjęciu Dotocz): X, Z 0,387C/s, «a „_ © 0,137C/5, 
t„ x [880-1150] s, h, = 261”C, t, = 3590 s, v, = 902?C, B=917,7'C, az 
= (,16”C/s, A = 721'C, a, — 0,22?C/s. 

Próbki początkowej fazy odpowiedzi skokowej — charakterystyki roz- 
grzewu oraz wyniki obliczeń n, podano w tabeli 6. Po odrzuceniu tylko 
jednego wyniku ny, < 2, uzyskuje się Średnio hy = 2,466, co oznacza n=3. 
Na podstawie zależności (41) można oszacować wzmocnienie k = 1440'C. 
Biorąc pod uwagę wyniki obliczeń T na podstawie związków (15) i (19) oraz 


uwzględniając [19] prawidłowość T  t„/(n— 1) = [440-575] s, można zało- 


Tabela 6. Dane pomiarowe i wyniki obliczeń nę pieca PM-6/1100 dla 50% P,, 


EEE EET 


2,671 | 2,463 | 2,656 
24,227 |2,0653 


O 0 JO U R W W 
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żyć T =518 s. Pozwala to na przyjęcie dla pieca pustego 50% mocy modelu 


1440 


8 (5) 55187 


Zakładając Tmin = 240 s, uzyskuje się A = 0,3 i model (8) pieca pustego 50% 
mocy | 
1440 


Mile LF 
(5) = 75240) (A 5800) (1 + s2660): 


co daje na podstawie równania (43) temperaturę końcową modelu Dymod = 
= 852 C. 

Powyższe modele dla pieca pustego zasilanego 50% mocy grzejnej mają 
charakter pomocniczy. Ważniejszy jest problem modelowania tego obiektu 
zasilanego 100% mocy znamionowej. W tabeli 7 podano próbki uzyskane 
z charakterystyki rozgrzewu po odjęciu temperatury otoczenia dla począt- 
kowej fazy tej charakterystyki pieca zasilanego pełną mocą znamionową oraz 
wyniki obliczeń no. Poniżej zestawiono też pozostałe dane możliwe do 
odczytania z charakterystyki rozgrzewu: «u, = 0,7'C/s, a„_ m 0,13”C/s, t, 
= [789--1030] s, h,„ = 358'C, t, = 1800 s, v, = 830,8'C, A = 649,6'C, 04 z 
= 0,52'C/s, B = 850,67C, az £ 0,42”C/s, hy = 2,726, stąd n=3. Parametry 
charakterystyki stygnięcia podane wcześniej są w tym przypadku oczywiście 
aktualne. 


Tabela 7. Dane pomiarowe i wyniki obliczeń ny pieca PM-6/1100 dla 100% EP,, 


BOO 


2,28 | 2,245 | 2,505 | 2,530 
2,081 | 2,762 | 2,773 


Na podstawie wzoru (41) oszacowane wzmocnienie wynosi k = 16957C. 
Obliczając średnią stałą czasową T wg założenia (15) i (19) oraz uwzględniając 
także zależność T = t„/(n— 1), po uśrednieniu przyjąć można T = 450 s. Zatem 
model 


; 
3 
4 
5 
6 
7 
8 
9 


1695 


OZ E5450* 


Logiczne jest przyjęcie T,;, = 240 s, tj. identycznej jak dla tego pieca zasilane- 
go 50% mocy znamionowej. Wtedy 4 = 0,42 i uzyskuje się model o różnych 
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stałych czasowych w postaci 


1695 
(1+ s240)(1 + s571)(1+ s1360): 


Daje to po zastosowaniu wzoru (43) temperaturę końcową modelu dla czasu 
t, = 1800 8, Dymoa 5 823'C, co niewiele różni się od v, rzeczywistej. 

Na przykładzie pieca PM-6/1100 można też zilustrować szacowanie chara- 
kterystyki statycznej, gdyż są w tym przypadku dostępne dane dla 50% i 100% 
mocy znamionowej oraz początki charakterystyk stygnięcia pieca od tem- 
peratur końcowych v,. Na podstawie zależności (40) można oszacować sto- 
pień nieliniowości u, ponieważ ii = 0,5 (50% mocy) odpowiada x = k/k = 
= 1440/1695 = 0,849, więc i = 0,177. Natomiast z pomiarów a, i a„_ na 
podstawie równości (35) uzyskuje się u = 0,183. Przyjmując i = 0,18, można 
narysować, opierając się na zależności (34), aproksymację charakterystyki 
statycznej pieca, jak to pokazano na rycinie 3. 


M (5) z 1 = 0,42. 


0,5 I 


Ryc. 3. Zidentyfikowana charakterystyka statyczna pieca PM-6/1100 
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W omawianym przykładzie zmierzono h, = 358”C. Zakładając, że h, = 
= (0,8k i stosując wzór (34) można obliczyć odpowiadające temu u dla 
nagrzewania i stygnięcia, a następnie na podstawie zależności (39) można 
wyznaczyć wartości pochodnych charakterystyki statycznej w tych punktach. 
Uzyskuje się odpowiednio: 

u, = 0,0476, > = 3,682  1u_ =0,426, dk = 0,651, 
0u 0u 
co daje na podstawie równania (38) u = 0,178. Należy zatem uważać, że w tym 
konkretnym przykładzie charakterystyka nieliniowa (34) dobrze modeluje 
rzeczywistą. Zauważmy, że do jej wyznaczenia niezbędna jest znajomość tylko 
stopnia nieliniowości j. Znając wartość k, można z tej charakterystyki przejść 
na konkretną wartość k dla dowolnych wartości procentowych mocy zasilają- 
cej piec w stanie ustalonym. Tu pod pojęciem stanu ustalonego rozumie się 
także wartości średnie sygnałów wejściowych i wyjściowych w czasie sterowa- 
nia piecem w przypadku regulacji dwupołożeniowej temperatury pieca lub 
wartości średnie sygnałów w przypadku regulacji ciągłej. W pracy [20] 
zaprezentowano oszacowanie modeli S(s) i M(s) omawianego w tym przy- 
kładzie pieca PM-6/1100 z wsadem. Należy zauważyć, że modele dynamiczne 
oszacowane prezentowaną metodą na podstawie danych pomiarowych stano- 
wią liniowy kompromis aproksymacji nieliniowego obiektu rzeczywistego dla 
konkretnej wartości mocy zasilającej ten obiekt elektrotermiczny w stanie 
ustalonym. Nie były one budowane jako składowe elementy struktury z ryci- 
ny 1. Wydaje się, że dobrym sposobem badania tak postawionego problemu 
byłyby symulacje, do których punkt wyjścia stanowiłyby oszacowania charak- 
terystyki statycznej 1 wyniki uzyskanych aproksymacji liniowych dynamiki. 

5. Prostyobiektlaboratoryjny w postaci mosiężnego pręta Z) 14 mm, 
długości 300 mm grzanego z jednego końca. W odległości 100 mm od dru- 
giego końca umieszczono w otworze pasowanym termoelement płaszczowy 
o Z I mm. Obiekt ten pozwala na zdjęcie pełnej odpowiedzi skokowej od tem- 
peratury otoczenia do temperatury maksymalnej, ustalonej. Temperatura była 
rejestrowana co 10 s, a czas rejestracji odpowiedzi wynosił 70 min. Zarejest- 
rowano też pierwsze 30 min stygnięcia. Tym sposobem uzyskano możliwość 
wyznaczenia z wykresu odpowiedzi parametrów modelu K(s) w postaci: 


5:00" ER 
IPS225" 


Wzmocnienie podano tu w mV, a opóźnienie zastępcze T, i zastępczą stałą 
czasową T, w minutach. Ponadto określono t,„ = 8,6 min dla grzania i t„ = 
= 9,5 min dla stygnięcia. W tabeli 8 podano wartości próbek początkowej fazy 
zarejestrowanej odpowiedzi skokowej oraz wyniki obliczeń. Dla czasu od O do 
140 s wartość odpowiedzi była zerowa. W tabeli, poza pierwszym fragmentem, 
podano próbki co 60 s oraz przyjęto, że opóźnienie, które tu wyraźnie 
występuje, jest równe 150 s = 2,5 min. 


K (s) z 
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Tabela 8. Dane pomiarowe oraz oszacowania parametrów obiektu laboratoryjnego 


Zidentyfikowane parametry 


1,49 1,59 1,66 | 155 236 310 383 446|4,76 5,28 5,64 5,89 
1,76 588 698 761 2,38 2,54 2,66 
743 882 918 2,07 2,16 
1075 1029 1,90 
991 


Odrzucając wartości T większe od 900 s uzyskuje się średnio dla grzania 
T= 516 s = 8,6 min. Podobnie, odrzucając w tym przypadku wszystkie warto- 
Ści pierwszego rzędu w tabeli, uzyskuje się średnio b, = 2,2:10 7”, co daje 
średnie wzmocnienie (z uwzględnieniem T = 8,6 min) k = 5,8 mV. Ponieważ 
wartość no jest bliska liczbie 2, należy się więc spodziewać niewielkiego 
rozrzutu (A £ 1). Przemawia także za tym równość t„=T [19]. 

Na podstawie końcowej fazy zarejestrowanej odpowiedzi można oszacować 
Imax. W publikacjach [18, 19] podano zależność przybliżoną: 


At 
k k=h(t"" 
k—h(t, + At) 


Tax m 


Zarejestrowana pełna odpowiedź skokowa umożliwia określenie: k = 
= 5,88 mV, h(t) = 5,79 mV, h(t, + At) = 5,85 mV, At = 10 min, co podstawia- 
jąc do powyższego wzoru pozwala oszacować T,,, % 9,1 min. Stosując układ 
(20) szacuje się A = 0,898. Wynik ten potwierdza przypuszczenie o bardzo 
małym rozrzucie stałych czasowych. Stąd też różnica pomiędzy wzmocnieniem 
k a identyfikowanym k będzie też niewielka. Na podstawie wyrażenia (17), dla 
n=2 


-(1+4) 
+ zU10 


3 =8,87.7 INV, 
n 
Wynik ten zgadza się prawie idealnie z wartością wzmocnienia k odczytaną 
z zarejestrowanej odpowiedzi h(t). Zatem przybliżony model liniowy pręta ma 
postać: 


Ę $.877e 7%? 


ARN SRG? 


Różni się od modeli przedstawionych wcześniej w przykładach tym, że 
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wystąpiło tu czyste opóźnienie, które zidentyfikowano na podstawie od- 
powiedzi. Można teraz zapytać, czy na podstawie ww. transmitancji, zastę- 
pując jej część bez opóźnienia 5,877(1+s8,6) * modelem K(s), uzyskamy 
zgodność tak otrzymanego modelu zastępczego z modelem zidentyfikowanym 
na początku tego przykładu z pełnej odpowiedzi. Stosując odpowiednio 
zależności podane w pracy [19] uzyskuje się 1, = T,, = 21,56 min, t,, z 
= 2,9 min, co pozwala oszacować T, = T+T,, = 5,4 min. Wyniki te są bardzo 
bliskie wyznaczonym na podstawie pełnej odpowiedzi skokowej. 

Początek zarejestrowanej charakterystyki stygnięcia umożliwia oszacowa- 
nie zastępczej stałej czasowej dla stygnięcia równej I,_ = 33,5 min, a to z kolei 
pozwala oszacować nieliniową charakterystykę statyczną dla określonego 
stopnia nieliniowości 


jA 
u=— 0,67. 


4. ZAKOŃCZENIE 


Należy podkreślić, że omawiana metoda umożliwia przybliżoną identyfika- 
cję modeli procesów, w szczególności elektrotermicznych, jako obiektów 
sterowania i regulacji. Przytoczone przykłady potwierdzają jej efektywność. 
Ważniejsze wnioski można streścić w kilku punktach. 

Przyjmując, że początek odpowiedzi skokowej to jej fragment od t = 0 do 
t = t„, podstawową zaletą metody jest prostota i szybkość uzyskania wyników 
identyfikacji, już bowiem w momencie t xt, znany jest rząd i pozostałe 
parametry modelu S (s). Rząd jest szacowany z dużą pewnością i jeżeli przyjąć 
operowanie rzędem całkowitym, bezbłędnie przy znajomości opóźnienia. 
Pozostałe parametry identyfikowane są z mniejszą pewnością — mogą 
wystąpić duże błędy. Poprawę można uzyskać dzięki konfrontacji wyników, 
ponieważ oszacowanie T i k jest możliwe przy znajomości dodatkowo t, na 
3 sposoby. Powyższe potwierdzają badania praktyczne. 

Szczególna przydatność praktyczna prezentowanej metody dla identyfikacji 
procesów elektrotermicznych wynika ze specyfiki rezystancyjnych pieców 
elektrycznych oraz z ogólnie przyjętego sposobu ich badania. Charakterystyka 
rozgrzewu pieca [6], która jest w rzeczywistości dużym fragmentem od- 
powiedzi skokowej, może być podstawą oszacowań dynamiki pieca, a fakt jej 
zakończenia przy wartości temperatury znamionowej pieca może być bardzo 
pomocny w oszacowaniu wzmocnienia k oraz maksymalnej stałej czasowej 
Imax. Ponadto, ponieważ zwykle znana jest dynamika czujnika temperatury 
stanowiącego integralną część obiektu, można znać wartość 7,,,,. Umożliwia to 
oszacowanie rozrzutu A i przejścia od modelu S(s) do modelu M (5), który jest 
modelem dokładniejszym. 

Przydatność tej metody w elektrotermii wynika także z powodu niemoż- 
liwego praktycznie bezpośredniego pomiaru wzmocnienia k obiektów elektro- 
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termicznych. Ponadto, bardzo powolne obiekty elektrotermiczne przy za- 
stosowaniu metod statystycznych identyfikacjj wymagają bardzo długiego 
czasu obserwacji. Prezentowana metoda daje wyniki już najpóźniej w czasie t,. 

Oszacowanie czasu t„, w którym charakterystyka skokowa osiąga punkt 
przegięcia, jest możliwe technicznie metodą próbkowania odpowiedzi skoko- 
wej. Znając t„; można sprawdzić oszacowanie 4 lub obliczyć A numerycznie 
[19]. Można także znając t, skonfrontować oszacowania T oraz k. 

Wpływ nieliniowości statycznej na wyniki oszacowania parametrów modeli 
nie jest znaczący z uwagi na to, że podstawą oszacowań są pary próbek 
początku odpowiedzi, położone ponadto blisko siebie. 

Wpływ zakłóceń losowych, nałożonych na sygnał użyteczny, jest teoretycz- 
nie duży z uwagi na możliwy niekorzystny stosunek szumu do sygnału 
użytecznego. Może ten wpływ być znacznie zmniejszony przez uśrednienia 
wyników obliczeń w serii, po odrzuceniu oszacowań skrajnych. Można też 
zastosować metodę filtracji cyfrowej [22] lub analogowej. 

Możliwe jest proste bezpośrednie oszacowanie przybliżonego modelu 
dyskretnego dla danego okresu próbkowania w czasie sterowania, co dodat- 
kowo pozwala na niezależne oszacowanie wzmocnienia, umożliwiając ewen- 
tualną konfrontację wyników. Tak uzyskany model dyskretny może być 
przydatny w sterowaniu cyfrowym. 

Bardzo ważne dla oszacowania stopnia nieliniowości jest to, aby skoro już 
zdjęto charakterystykę rozgrzewu, zarejestrować także początkowy fragment 
stygnięcia. 
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STANISŁAW SKOCZOWSKI 


IDENTIFICATION AND MODELLING OF ELECTROHEAT 
PLANTS ON THE BASIS OF INITIAL PHASE OF HEATING-UP 
CHARACTERISTIC 


Summary 


The theoretical basics and experimental research results of a new deterministic method of 
simplified models used in electric resistance furnaces are presented in the paper. Assuming linearity, 
the method allows to identify models on the basis of discrete samples of the step response's initial 
phase, which is synonymous with a start-up period of furnace heating. 

Obtained very quickly identification results can be used, e.g. in auto-tuning of a furnace's 
temperature controls. 
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CTAHHCJIAB CKOUOBCKIA 


MIEHTHOKKANMA H MOJNEJJMPOBAHME 
DJIEKTPOTEPMMUECKHX OBBEKTOB HA OCHOBE 
HAUAJIBHOŃ QA3bIl XAPAKTEPACTHKI IIOJIOTPEBA 


Pe3roMe 


B pa6oTe oroBapuBaroTca TeopeTuueckHe OCHOBbI HM DpE3YJIbTATbI JKCIIEDUMEHTAJIBHBIX 


UCCJIETOBAHHA HOBOTO METOJA HNeHTUDHKANKM YIIPpOLIEHHbBIX MOJIEJEŃ HPpH IpUMEHeHHH NIA 
pPe3ACTAHIIAOHHBIX Ieueń. IIpA NpeNNoJloxeHiM JIMHEĄHOCTH METOJĄA, UTO IIO3BOJIACT HUN1EH- 
TUQUUHMPOBATb MOJEJI5 Ha OCHOBE /MHCKPETHBIX ACHBITAHHŃ, HaUAJIBHOŃ (Da3bl CKAUKOOÓPAZHOTO 
OTBETA, OJIHO3HAUHOTO C XAapaKTEpHCTHKOŃ IOJĄOFPEBA IIeUH. 


IlojryueHHbie TakAM OÓpa30M MTHOBeHHbie DpE3YJIbTATbI AUNCHTHQHKANAA MOTYT ÓbiTb 


ACIIOJIB3OBAHbI, HAHDHMEDp, NJIA CAMOHACTpOŃKH peryJiaTopa TeMliepaTypkl IIeYH. 
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Streszczenie: Przedstawiono uproszczoną analityczną metodę określania 
pola temperatur w pasmowych strukturach cienkościennych. Metoda ta oparta jest 
na założeniu jednowymiarowego przewodzenia ciepła w elementach struktury. 
Omówiono sposoby rozwiązania zagadnienia dla dwóch wariantów: w wariancie 
A pojemność cieplna czynników omywających elementy struktur była nieskończenie 
wielka, w wariancie B pojemność cieplna strumieni czynników omywających 
elementy struktur miała skończoną wartość. Dla wybranej struktury dokonano 
porównania rozwiązania uproszczonego jednowymiarowego z rozwiązaniem nie- 
uproszczonym dwuwymiarowym, na tej podstawie podano zakres stosowania 
rozwiązań jednowymiarowych. 


Słowa kluczowe: równanie przewodnictwa — warunki brzegowe — ustalone 
przewodzenie ciepła — jedno- i dwuwymiarowe pole temperatur. 


1. WPROWADZENIE 


Przy ogrzewaniu różnego rodzaju płyt, den i ścian zbiorników stosowany jest 
sposób prowadzenia czynnika grzewczego w profilach półcylindrycznych lub 
kątowych przyspawanych do jednej z powierzchni podgrzewanych elementów. 
Ułatwia to czyszczenie zbiorników, jak również pozwala w narzucony techno- 
logicznie sposób kształtować drugą powierzchnię płyty. Na podobnej zasadzie 
odbywa się chłodzenie skraplaczy lodówek. Czynnik używany do podgrzewania 
(chłodzenia) może zmieniać stan skupienia (czynnik o nieskończonej pojem- 
ności cieplnej strumienia cieplnego) lub może zmieniać się jego temperatura 
(w przypadku skończonej pojemności cieplnej strumienia czynnika). 

W wymienionych wyżej rozwiązaniach ścisłe określenie ilości przekazywane- 
go ciepła jest na ogół bardzo złożone, a w niektórych przypadkach wręcz 
niemożliwe. Z tego powodu przy konstrukcjach złożonych, wykonywanych 
z powtarzalnych cienkościennych elementów, można zastosować przybliżone 
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metody, uwzględniające tylko jednokierunkowe przewodzenie ciepła [1, 6]. 
W niniejszym opracowaniu przedstawiono przybliżoną metodę określania pola 
temperatur w strukturach pasmowych. Struktury pasmowe to takie struktury 
cienkościenne, w których wydzielić można elementy w kształcie pasm, różniące 
się warunkami wymiany ciepła. Przykład takiej struktury przedstawiono na 
rycinie 1. 


N CZYNNIK GRZEWCZY 


Szerokość _ płyty 


Ryc. 1. Przykład struktury pasmowej 


Dla jednej wybranej struktury dokonano porównania rozwiązania uprosz- 
czonego jednowymiarowego z rozwiązaniem dokładnym dwuwymiarowym. 
Na tej podstawie określono zakres stosowania metody uproszczonej. 


2. DWUWYMIAROWE POLE TEMPERATURY 
W PASMOWYCH ELEMENTACH CIENKOŚCIENNYCH 


2.1. SFORMUŁOWANIE ZAGADNIENIA 


W strukturach pasmowych można wydzielić charakterystyczny element, 
taki, który z jednej strony częściowo jest nagrzewany, z drugiej zaś oddaje 
ciepło. Schemat takiego elementu przedstawiono na rycinie 2, a tok obliczeń 
podano poniżej. 


Bi,0=0 


Ryc. 2. Schemat charakterystycznego elementu 
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Po wprowadzeniu wielkości bezwymiarowych: 


W=Ż możę Bis; (2.1) 
ć c A 


d T—t 
Ąd== D=—, OR LEE 


c c z +b6 


(2.2) 


warunki brzegowe dla omawianego zagadnienia można opisać następująco: 


60(X,Z)| 
00.2] _, 2) 


00 (X, Z) 


= (0, ? 
2 la. (2.4) 
00 (X, Z) | 
O zzRBiG. 5 
„zag b I (2.5) 


Dla Z = D można sformułować trzy wersje warunku brzegowego, w zależ- 
ności od określenia wielkości charakteryzujących wymianę ciepła. Można 
przyjąć, że wnikanie ciepła następuje przez część powierzchni określoną 
w przedziale od 0 do 4. Pozostała część powierzchni w przedziale od A do 
C=| jest adiatermiczna (560/02 = 0). 


2.1.1. WARUNEK BRZEGOWY PIERWSZEGO RODZAJU 


Wymiana ciepła z czynnikiem w przedziale od 0 do A odbywa się przy 
określonej temperaturze powierzchniowej T = T,(© =0), stąd 


00 (X, Z 
O (X, ZJr-pk(X + i [1-k(X)] = k(X), (2.6) 

gdzie i 
k(X0lxE86=1, k(XJ|z=A = 0. (2.7) 


2.1.2. WARUNEK BRZEGOWY DRUGIEGO RODZAJU 


Doprowadzenie ciepła z ośrodka w przedziale od 0 do A odbywa się przy 
zachowaniu stałej gęstości strumienia cieplnego: 


00 (X, Z) 


1 
Ż = -k(X)q. (2.8) 


z=b A 


2.1.3. WARUNEK BRZEGOWY TRZECIEGO RODZAJU 


Wymiana ciepła z czynnikiem w przedziale od 0 do A odbywa się według 
prawa Newtona, dlatego 
00 (X, Z) 
OZ 
W zależnościach (2.8) i (2.9) funkcja k(X) określona jest równaniem (2.7). 


= k(X)Bi;[1—0(X, Y)|z=pl. (2.9) 


Z=D 
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2.2. ROZWIĄZANIE 


Równanie przewodnictwa cieplnego dla stanu ustalonego, bez wewnętrz- 
nych źródeł ciepła, sprowadza się do równania Laplace a: 


26 00 


jyztzzz = 0. (2.10) 


Dla warunków brzegowych, określonych równaniami (2.3), (2.4) i (2.5) roz- 
wiązanie przybiera postać: 


O(X, Y)= 27 „Ccos(znX) JZĘENCZJ (2.11) 


Dla Z=D rozwiązanie przybiera różną postać, w zależności od warunku 
brzegowego. 


2.2.1. DLA WARUNKU BRZEGOWEGO PIERWSZEGO RODZAJU 


Z warunku brzegowego określonego zależnością (2.6), po uwzględnieniu 
równania (2.11) i dalszych przekształceniach, otrzymuje się: 


f cosh (znD) [k (X)—Bi(1—k(X))]+ 
). E,cos(nnX) 4 Bi =k(X). (2.12) 
n=1 |IZZZZJ 


Jeżeli, po uporządkowaniu, wyrażenie w klamrze opiszemy przez 


p„(x) =cosh (znD) [k (X) —Bi(1 —k(X))] — 
-siah (mD) | KC) rn(1-ktv) | (2.13) 
to równanie przybiera postać: 
> E,cos(anX) p, (x) = k(X). (2.14) 
Po obustronnym wymnożeniu zależności (2.14) przez cos(xnX) oraz scał- 


kowaniu otrzymano układ nieskończonej liczby równań, umożliwiający wy- 
znaczenie stałych E,: 


d1,1 d1.2 d1,—1 dn E, bi 

da1 da, 2 dajn-1 d2 jn E; bą 
PREY 7 A ER WERE FCC X zj = w "(Żi5) 
dn= 11  dn-1,2 dn-1n-1 dn-1mn Bza Ba 

a, 1 a, 2 dnn-1 dp.n n b, 
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gdzie 
1 
a,, = | cos(znX)cos (xkX)p,(X)dX. (2.16) 
0 


Elementy macierzy współczynników określone są wzorem: 


_ Jsn[z4(n—k)] sin[xrA(n+k)] 
nk | 2(n=k)n 2(n+k)n 


x |oosh (znD) (1 + Bi) — sinh (znD) Ę + m || (2.17) 


Wyrazy b, określa się z zależności: 


b, = —(k(X) cos (nkX)dX, (2.18) 


O 


przy czym po dalszych przekształceniach 
l 
b, = ——sin(xkA). (2.19) 
nk 


Po rozwiązaniu układu równań (2.15) i określeniu współczynników pole 
temperatury określa się ze wzoru (2.11). 


2.2.2. DLA WARUNKU BRZEGOWEGO DRUGIEGO RODZAJU 


Uwzględniając warunek brzegowy (2.8) w równaniu (2.11), można je prze- 
kształcić do postaci: 


ŻE „cos(nnX) sh (anD)—>cosh tny) F -k(X)?. (2.20) 


Po wykorzystaniu ortogonalności funkcji własnej i dalszych przekształceniach 


2q sin(unA) 


E = (2.21) 


n 


1(an)? sim (unD) — cosh (np) | 


Pole temperatury określa się z zależności (2.11) po wprowadzeniu doń 
wyrażenia (2.21). 


2.2.3. DLA WARUNKU BRZEGOWEGO TRZECIEGO RODZAJU 


Po uwzględnieniu warunku brzegowego (2.6) w równaniu (2.11) można je 
przekształcić do postaci: 


DE „cos(unX)p,(X) = —k(X)Bi,, (2.22) 


3 — Szczecińskie Roczniki Naukowe t. XI, z. 1 
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gdzie 
B; 
p,(X) = sinh (unD) |m + k(X) Bi, 4 —cosh (anD) [Bi+k(x)Bi,]. (2.23) 


Po obustronnym wymnożeniu zależności (2.22) przez cos(anX) oraz scał- 
kowaniu otrzymuje się również układ nieskończonej liczby równań, umoż- 
liwiający wyznaczenie stałych E„. Zagadnienie to sprowadza się do zagadnienia 
z warunkiem brzegowym pierwszego rodzaju. Elementy macierzy współczyn- 
ników definiowane są wzorem: 


s . |sn[zrA(n—k)] sn[xA(n+k)]|| . Bi 
ax = —Bi, 2m=%a — TY: si (znD) z* + cosh (Dy) 


(2.24) 


Wyrazy b, określa się z zależności: 
B; 
b, = — —Ysin (zkA). (2.25) 


Po rozwiązaniu układu równań (2.14) pole temperatury ustala się na podstawie 
zależności (2.11). 


3. ELEMENTY CIENKOŚCIENNYCH STRUKTUR PASMOWYCH 


3.1. ZAŁOŻENIA 


Zakłada się, że elementy struktury są na tyle cienkie, że można pominąć 
w nich rozkład temperatury w przekroju poprzecznym. Analizując struktury 
pasmowe, można wyróżnić dwa podstawowe typy elementów: 

typ I — pasy, dla których na obu brzegach zachodzi przekazywanie ciepła 
z sąsiadującymi elementami (temperatura określona na obu brzegach), 

typ II — pasy, których jeden brzeg jest adiatermiczny, przekazywanie ciepła 
zachodzi na jednej stronie (temperatura określona z jednej strony, druga strona 
adiatermiczna). 

Dla pasów tych można poczynić następujące założenia: stała wartość 
współczynnika przewodności cieplnej, stała wartość współczynnika przej- 
mowania ciepła, jednowymiarowe pole temperatury, stałe temperatury otocze- 
nia i na brzegach pasów, rozkład pola temperatury ustalony w czasie. 

Przy tychże założeniach strukturę pasmową można podzielić na elementy 
pokazane na rycinie 3. Przyjmuje się dla każdego elementu ij lokalny 
niezależny układ współrzędnych x,. Dla sformułowanego zagadnienia brzego- 
wego określenie pola temperatury sprowadza się do rozwiązania równania 
różniczkowego [3]: 

2 T.. xf OL.:; 
R = go CROFT 64) 
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Ryc. 3. Schemat elementu struktury pasmowej 


Równanie powyższe można przedstawić w postaci: 


d* T,(%%) 

DE A - — Ly Ty) +9; = 0, (3.2) 
XK 
gdzie 
AQij F Ugij ; 

g TORA Świ: 3.3 
kij d,, U ( ) 
ojj loij F Zgij i gij. 3,4 
di, — dzaj ( * ) 


3.2. ROZWIĄZANIE DLA ELEMENTU TYPU I 


Na brzegach elementu typu I spełnione są następujące warunki brzegowe: 


TEG =0 r” H, (3.5) 
BG|<=1, 75: (3.6) 


Po podstawieniu warunków brzegowych do rozwiązania ogólnego otrzy- 
muje się 


Jij 
T;,(%%) zę ( U) cosh (Xx L;) + 
ij 
Jij Ji 
( j-u)-(q- du cos (1, NIM) R 
— ———-—=sinh (x 4/ 1) +. (3.7) 
sinh (l,, 4/ ,,) Mij 
Znając rozkład temperatur w elemencie, można określić strumienie cieplne na 


brzegach pasów. Strumień cieplny przekazywany do brzegu o temperaturze T; 
(dla x, = 0) określony jest zależnością: 


QLij F UAT T;+D;; „AJ T;+4;; (3.8) 


=bddyyj Hgetelilso/ Hs). (3.9) 


gdzie 
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1 — Ewa hy, Pij 

7 sinh (1, / u) 

bd,,A,;g,, [cosh ( (1, /u,)— 3.11) 


di,i 
s" NH sinh (l, R 


Strumień cieplny przekazywany do brzegu o temperaturze T,(x, = l,,) okreś- 
lony jest zależnością: 


(3.10) 


Qyji = Ai DD jj TQ ji (3.12) 


Aj = —bdjj A, / u, ctgh (l, 4/ u), (3.13) 
b, = —bd,; 4, +/p,, [sinh ( (l; / 14) cosh (1, /p4;) ctgh ( l; /45)]. (3.14) 


43% _ bdzj kuj g,j [cosh (1, / Hs) 1] NHs)—1] (3.15) 
eg SZ sinh (1, / u.) 


3.3. ROZWIĄZANIE DLA ELEMENTU TYPU II 


gdzie 


Na brzegach elementu typu II spełnione są następujące warunki brzegowe: 


T;;(X%)|x4=0 = b, (3.16) 
T:. 
ATC) = (3.17) 
dy, Xki = lij 
Rozwiązanie ogólne przyjmuje postać: 
(1- Zu cos [(1 xy | g;, 
T,(24) = — +—. (3.18) 
ach cosh (1,, 4/1) Hij 
Strumień cieplny określony jest zależnością: 
Qi; = A; jj Ti + Qi, ij? (3.19) 


gdzie 


—bdyj A, / 4 tgh (li / kj). (3.20) 
bd,,4,,g; 
di,ij — ma gh (l;, i kt). (3.21) 


ij 
3.4. PODSUMOWANIE 


Dla elementu między dwoma węzłami o dowolnych numerach i,j, na 
występujący w zależnościach na strumień cieplny na brzegach wyraz wolny, 
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zachodzi równość: 


diij = djij 7 dij: (3.22) 
Wtedy strumienie cieplne można określić ze wzorów: 
Qi; 7 dy TG +b; Ti+4j> (3.23) 


jeżeli strumień w elemencie ij skierowany jest ku brzegowi o temperaturze 
1,, oraz 


Qi = Gi +0 GŁ dy (3.24) 


*J 


jeżeli strumień w elemencie ij skierowany jest ku brzegowi o temperaturze 7,. 
Dla elementów utwierdzonych w jednym końcu w wyrażeniu na strumień 
cieplny nie występuje temperatura końca nie utwierdzonego. 


4. POLE TEMPERATURY PRZY PODGRZEWANIU PASMA 
CZYNNIKIEM O NIESKOŃCZONEJ POJEMNOŚCI CIEPLNEJ 
(WARIANT A) 


4.1. PODSTAWY METODY 


W rozważanych strukturach wydzielone pasy łączą się ze sobą w liniach 
węzłowych zwanych w dalszych rozważaniach „węzłami”. Równanie bilansu 
energii dla poszczególnych węzłów można zapisać w postaci: 


2 Qi; w Q;, (4.1) 


gdzie Q, oznacza strumień doprowadzany do i-tego węzła zewnętrznego. Dla 
węzłów wewnętrznych Q,=0. 

Z równania bilansu (4.1) dla poszczególnych węzłów, wykorzystując dla 
każdego pasa łączącego się w węźle zależności na strumień ciepła (3.23) i (3.24), 
otrzymuje się układ równań algebraicznych, w którym niewiadome są tem- 
peratury. Układ ten umożliwia obliczenie temperatur w poszczególnych 
węzłach. A znając ich wartości określić można rozkład temperatury w po- 
szczególnych pasach. 

W układach bardziej złożonych, składających się z dużej liczby węzłów 
1 elementów, dla ułatwienia obliczeń z wykorzystaniem elektronicznej techniki, 
układ równań należy sprowadzić do ogólnej postaci macierzowej [3]. W tym 
celu sprecyzowany zostanie bliżej formalizm opisu takiego układu. 

Do dalszych rozważań przyjmuje się układ składający się z n węzłów, 
w których istnieje oddziaływanie cieplne pomiędzy wszystkimi węzłami, 
łączącymi się pasami. Przy takim połączeniu n węzłów (każdy z każdym) 
w układzie będzie występowała maksymalna liczba n!/[2(n—2)!] elementów. 
Do każdego węzła z n węzłów, które traktowane są jako węzły zewnętrzne, 
doprowadzany jest strumień Q,. Na każdy element mogą oddziaływać dwa 
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ośrodki o różnych temperaturach t g,, oraz t, ,,. Jeżeli między dwoma węzłami 
występuje więcej elementów, to do obliczeń można wprowadzić element 
zastępczy. Przykładowo, jeżeli między dwoma węzłami występują dwa ele- 
menty, to całkowity strumień ciepła w elemencie zastępczym ij, skierowany ku 
węzłowi o temperaturze T;, 


Qży = Qi; + Qi; (4.2) 


natomiast całkowity strumień ciepła w elemencie zastępczym ij, skierowany ku 
węzłowi o temperaturze T5;, 


Qfyi = Qi + Diji- (4.3) 


Strumienie ciepła występujące w tych zależnościach określa się na pod- 
stawie podanych wcześniej wzorów. 

Dla układu składającego się z n węzłów, wykorzystując równania bilansu 
dla poszczególnych węzłów oraz odpowiednie wyrażenia na strumienie cieplne 
skierowane do poszczególnych węzłów, otrzymuje się układ n równań alge- 
braicznych liniowych. Układ tych równań można zapisać w postaci macie- 
rzowej: 


C1,1 UE Bei Din Ti d, 
da1 C2,2 bajn-1 Dan T, d, 
MOE ERP MGLE BT. PU I R 3 | EEE (4.4) 
Adn—1,1  dn-1,2 EA NLREZ A f pode dzi 
An, Xq,2 Ann 1 Cn n | d, J 
lub 
AxT=D, (4.5) 
gdzie 
A — macierz współczynników równania, 
T — wektor kolumnowy niewiadomych, 


D — wektor kolumnowy wyrazów wolnych. 
Współczynniki występujące w macierzy A są określone przez zależności: 


Ci; = 4; +0D,; = » one Re, Ma 2, 22, Że, FE (4.6) 
j=f+1 j=ń 
yy „Ala J= 1, 2,34...,0—1 
oraz i=jJ+l,]+2,...,N (4.7) 
b,,=b,,, dla i=1,2,3,...,n—l 
oraz j=i+1,i+2,...,h. (4.8) 


Współczynniki w wektorze kolumnowym D określa się natomiast z za- 
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leżności: 


d=dr= Yeqy sdlai=1, ydyh. (4.9) 


j+1 
4.2. ZASTOSOWANIE METODY DLA STRUKTUR PASMOWYCH 


Największa struktura pasmowa, w której węzły są ze sobą połączone, każdy 
z każdym, wynosi 4. Na rycinie 4 przedstawiono ideowy schemat takiej 
struktury. 


Ryc. 4. Schemat pełnego układu czterowęzłowego 


Dla tej struktury podano poniżej macierz współczynników A oraz wektor 
kolumnowy wyrazów wolnych D. 


A1,1;2 +41,1;3 HF 41,1,4 by 4:2 by 1:3 by 1:4 
As S2,1,2 2,23 +A2,2;4 + Da, 2;1 ba, 2;3 ba, 2:4 
A3,1;3 a3,2;3 3,34 + b3,133 +03, 3;4 bz 3:4 j 
d4,1,4 d4 2:4 a4,2:4 bą ay +Dą 42 Hoa 433 


(4.10) 


Q; —(91:2 + 1:3 +41:4) 
— » +d>.2 + d>. 

D- Q+—(41:2 + da:3 2:4) | (411) 
Q> —(91:3 + qa.3 + q3:4) 


Q4 p (q4:4 + q2:4 nić q3:4) 


W strukturach pasmowych bardzo często występują układy szeregowe 
(w węźle łączą się tylko dwa pasma). W przypadku układu szeregowego 
n-węzłowego wektor kolumnowy niewiadomych T oraz wektor kolumnowy 
wyrazów wolnych określone są jak w równaniu (4.4), natomiast macierz 
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A współczynników równania przybiera postać: 


Paa o Ehe: 0 * 0 0 0 
dz, Ad27+b22 Das ... 0 0 0 
0 Az | dzą +D33 ser 0 0 0 
di MACOS MBP FSA AAA (4.12) 
0 0 0 Sj dą ży =25DRS245-2 oai 0 
0 0 0 a. An-1,n-2 Ara chB <w BR 
0 0 0 wsz 0 Gas4 b, 


Współczynniki macierzy 4, zależność (4.12), i wektora kolumnowego D można 
określić z następujących związków: 


o — EE: dl u ls:24 SLTAA —|] 
MOJETTZS 4 4 a: (4.13) 
oraz j=i+l, 
ma == SER dl = 2% Each 
by, = Dj ń : c ś (4.14) 
oraz j=i—l, 
d, = Q;—(4q;;-1 Fdii+ 1) dla l = l, Ż, 3 .... n. (4.15) 


Jeżeli węzeł początkowy i końcowy są węzłami adiatermicznymi, to układ 
n równań obejmuje niewiadome temperatury we wszystkich węzłach, z wyjąt- 
kiem węzłów adiatermicznych, których numeracja zaczyna się od cyfry zero. 

Na podstawie przedstawionych zależności można określić wartości tem- 
peratury w węzłach i rozkłady tewmperatury w pasach dla dowolnej struktury 
pasmowej. Poniżej przedstawiono rozwiązanie dla najprostszego układu, 
składającego się z dwóch pasów. Jest to układ szeregowy jednowęzłowy. 
Schemat takiego układu, odpowiadającego charakterystycznemu elementowi, 
przedstawiono na 'rycinie 5. 


olg | tg ol—0 Ag | tg 


Ryc. 5. Schemat układu szeregowego jednowęzłowego 


Równania dla tego układu przybierają postać: 
Q1.,1:0,1 = dyo1 i +44:01: (4.16) 
Q1.1:02 = 41:02 i €41:02> (4.17) 
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stąd 
T a __du;01 d102: (4.18) 
A1.01 HF A1.02 


5 PORÓWNANIE METODY DWUWYMIAROWEJ 
Z JEDNOWYMIAROWĄ DLA WARIANTU A 


"W układach szeregowych zbudowanych z powtarzalnych sekcji, w których 
temperatury czynników nie zmieniają się, a parametry odpowiadających sobie 
elementów są takie same, rozkłady temperatur w poszczególnych sekcjach 
powtarzają się. Schematyczny przekrój przez układ przedstawiono na rycinie 6. 


© oloor to | QT) or tor © 
WAJALGG666, 


Olgo1  too1 Olgo1=0 tgo1 


RR"MF PORE: SWR 


Ryc. 6. Schematyczny przekrój przez układ 


Powtarzalna sekcja odpowiada elementowi, dla którego określono dwu- 
wymiarowe pole temperatury. Z analizowanych warunków brzegowych naj- 
bardziej zbliżone do rzeczywistych są przedstawione w wariancie A. Dla tego 
wariantu dokonano obliczeń rozkładu temperatur dla kilku kompletów 
danych wejściowych. Przykładowy rozkład temperatury dla przypadku skraj- 
nego (skraplanie pary wodnej z jednej strony, z drugiej zaś przepływ powietrza 
z prędkością 14 m/s) przedstawiono na rycinie 7. Można zauważyć, że 
w przekroju poprzecznym temperatura zmienia się w niewielkim zakresie. 

W dalszej kolejności przeprowadzono obliczenia dla tych samych kom- 
pletów danych przy założeniu jednowymiarowego pola temperatury. Analizo- 
wany przypadek odpowiada najprostszemu układowi szeregowemu z jednym 
węzłem. 

W tabeli 1 zamieszczono obliczone dwiema metodami wartości temperatur 
w charakterystycznych punktach sekcji dla kompletu danych odpowiadające- 
mu przypadkowi przedstawionemu na rycinie 7. Maksymalna względna 
odchyłka wyniku dla rozwiązania jednowymiarowego w stosunku do wyniku 
rozwiązania dwuwymiarowego, dla wybranych charakterystycznych wartości 
x, nie przekracza + 1,96% dla x =0, + 2,86% dla x = 0,042 m oraz 3,66% dla 
x=0,175 m. 

W przypadku gdy liczba Biota liczona względem czynnika grzejnego 
Bi, < 0,5 i gdy stosunek liczb Biota względem czynnika grzejnego i grzanego 
Bi,/Bi > 20 bądź Bi < 0,01, maksymalne różnice między wartościami tem- 
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POLE TEMPERATURY 
zm] 


d=0,01 1 1029C__ 10090 957C 90?C 80% 70% 60% 509C 40?C 30?C 


1039C 


[o 


0,0 0,05 0,10 0,15 c=0,175m 


Ryc. 7. Pole temperatury i temperatura Średnia dla elementu płyty przy ogrzewaniu parą na- 
syconą dla c = 0,075 m, a = 0,042 m, d=0,01 m, ty = —40”C, £,=110'C, a = 36 W/(m* :K), 
= 1500 W/(m?*:K), 4 = 40 W/(m:K) 


peratury określonymi na powierzchni górnej i dolnej oraz wartością określoną 


z modelu jednowymiarowego względem temperatury odniesienia mogą wyno- 
sić około 20%. 


Tabela 1. Wartości temperatur w charakterystycznych punktach sekcji 


Współrzędne [m] 00 | 0,042 


Rozwiązanie dwuwymiarowe 
z=001lm=d 


z=0 
wartość średnia 


Rozwiązanie jednowymiarowe 


6. POLE TEMPERATURY PRZY PODGRZEWANIU PASMA 
CZYNNIKIEM O SKOŃCZONEJ POJEMNOŚCI CIEPLNEJ 
(WARIANT B) 


6.1. ZAŁOŻENIA 
Przyjmuje się, że w strukturze występują pasma ogrzewane czynnikiem 


o skończonej pojemności cieplnej [4, 5]. Czynnik przepływając oddaje ciepło 
i jego temperatura obniża się. 
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Na podstawie analiz zjawiska wymiany ciepła przez pas i opisanych 
w piśmiennictwie przyjęto następujące założenia do modelu matematycznego: 

— wymiana ciepła między czynnikiem a płytą oraz między płytą a ośrod- 
kiem o stałej temperaturze odbywa się na drodze wnikania ciepła, 

— temperatura t w przekroju poprzecznym czynnika jest wielkością stałą, 

— temperatura T na grubości płyty jest stała, 

— pomija się przewodzenie ciepła w płycie i czynniku wzdłuż pasm, 
pojemność cieplna czynnika jest stała. 


6.2. METODYKA ROZWIĄZANIA ZAGADNIENIA 


Wycina się ze struktury poprzeczny pasek o elementarnej szerokości dy, jak 
przedstawiono na rycinie 8. W pasku takim, analogicznie jak w przypadku 
struktury o stałej temperaturze czynników, wyodrębnić można węzły i elemen- 
ty przewodzące ciepło. 


(bQL/LLkokok LLL kok ALLAL LL LL fok l olL LA kok l kokok A LM 4d 4LLĄ LLL TG LSLLL 
J << 3 A uj 
= 


MZ WE 


Ryc. 8. Schemat płyty ogrzewanej czynnikiem płynącym w kanałach 


Przy pominięciu rozkładu temperatury na grubości elementu można, dla 
elementarnego paska o szerokości dy, założyć jednowymiarowe pole tem- 
peratury w poszczególnych elementach. W rozpatrywanym pasku o szerokości 
dy, dla wyodrębnionych węzłów o współrzędnych y, równanie bilansu energii 
ma taką samą postać, jak dla rozważanej struktury o stałej temperaturze 
czynników: 


Ż:;(9) = Q;,(y). (6.1) 
j 


Otrzymuje się analogiczny układ równań, liczący tyle równań, ile wyodręb- 
niono węzłów. Dla stanu ustalonego strumienie ciepła dopływające do danego 
węzła j od elementu ij można określić z ogólnej przedstawionej już zależ- 
ności (3.8): 


QW) = Gzy 7,0) +bz p 709) Fi (6.2) 
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W tym przypadku temperatury i strumienie cieplne są funkcjami współ- 
rzędnej y. Współczynniki a, b, q zależą od parametrów wejściowych: współ- 
czynników wnikania ciepła, przewodności cieplnej materiału, wymiarów linio- 
wych, temperatur otoczenia elementu i temperatury czynników. Przyjmuje się, 
że temperatury czynników są nieznane. 

Rozwiązując układ równań (6.1), otrzymuje się temperatury w węzłach, 
które określone są w funkcji temperatur czynników t,(y). Zależności na 
temperatury w węzłach mają więc postać: 


T;()) = Ż Bi: ti Łój. (6.3) 


Litera m, występująca nad znakiem sumy, oznacza liczbę czynników wpływają- 
cych na temperaturę w danym węźle. Współczynniki $,, i ó, są sumą 
odpowiednich współczynników a, b, 4, i podobnie jak one zależą od paramet- 
rów wejściowych. 

Na podstawie wartości temperatur w węzłach można określić temperatury 
elementów wzdłuż rozpatrywanego paska, jak również temperaturę Średnią dla 
każdego z nich. Temperaturę średnią poszczególnych elementów określa się 
z zależności: 

Lji 
» | 
1,(y) = m | TW, xy) dxy. (6.4) 
Jt A 
Po scałkowaniu powyższą zależność można przedstawić jak wyrażenie (6.3) 
w postaci: 


T;, (Y) = ż Hak (9) + 0;,. (6.5) 


Występujące w tym wyrażeniu współczynniki zależą również od parametrów 
wejściowych. 

Mając określone temperatury średnie dla wszystkich elementów paska, 
o elementarnej szerokości dy, można napisać równanie bilansu cieplnego dla 
poszczególnych czynników: 


.... . a. a a m a a a ; .. . a am a a .. . . „a a. a w , (6.6) 
+ W,_4dt,_, = 0-1 ( > [t,-1(9)— T,.:;0)] Jiji dy 
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Litera k, występująca nad znakiem sumy, oznacza liczbę elementów 
stykających się z danym czynnikiem. Znak + zależy od kierunku przepływu 
czynnika. Z równań bilansu otrzymuje się układ równań różniczkowych 
liniowych pierwszego rzędu niejednorodnych: 


k 


dt, E3 
ESA +2 A, ,t,(0)=D 


k 
tz+ ż A> „t,(y) = D» 
NIP OEM FUN A ANA (6.7) 


RAR 2. A; t;(v) =D 


1= 


Powyższy układ równań rozwiązuje się dla warunków brzegowych, którymi są 
temperatury wlotowe czynników: 


U (Wls=y+, — l1p 
lą (Vly=y>, = lap 
HRT" (6.8) 


PRZE (Potato = tn—1)p 
Ln1 (3|»= x, = tut 


W wyniku rozwiązania układu równań różniczkowych otrzymuje się funkcje. 
Funkcje te opisują zmienność temperatur czynników wzdłuż pasm. Znając 
rozkład temperatur czynników, określić można z równania (6.3) zmienność 
temperatur wzdłuż linii węzłowych i dalej obliczyć temperatury w dowolnych 
punktach elementu. 


6.3. WYNIKI OBLICZEŃ 


Obliczenia rozkładu pola temperatur wykonano dla płyty podgrzewanej 
czynnikiem o skończonej pojemności cieplnej o temperaturze wlotowej 7, przy 
współczynniku wnikania ciepła «,, płynącym w kanałach znajdujących się w jej 
spodniej części [ 2, 3 |, zaizolowanej termicznie od spodu oraz chłodzonej od góry 
przez czynnik o stałej temperaturze ty, przy współczynniku przejmowania ciepła 
a. Założono, że parametry i rozstaw kanałów były na całej szerokości płyty takie 
same. Czynniki w kanałach miały taką samą pojemność cieplną, a ich przepływ 
odbywał się w przeciwprądzie. Schemat płyty pokazano na rycinie 8. 

W płycie wyodrębniono powtarzalny element. Z przedstawionych powy- 
żej ogólnych zależności otrzymano dla tego elementu układ dwóch równań 
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różniczkowych, rozwiązano ten układ i dalej określono temperaturę w wy- 
branych punktach płyty. Wyniki obliczeń dla fragmentu płyty zawartej między 
kanałami przedstawiono na rycinie 9. 


23,399C 


48,29?C 
12,74 


0,1 0,2 0,265 x[m] 

Ryc. 9. Izotermy dla fragmentu płyty zawartej pomiędzy profilami dla t,(0) = t,, , (b) = 333 K, 
to = 233 K, W, =W,,, = 168 W/K, a, = 36 W/(m*:K), a, = 2000 W/(m*:K), a, = 0 W/(m*:K), 
Asian 7 Ai-q j 540 W/(m:K) d,;,, =d,,-1=0,017 m, b=127m, I,_,,= 0,083 m, I 
= (0,265 m 


ii+1 


7. WNIOSKI 


Przedstawiono analityczną metodę określania pola temperatury w pas- 
mowych strukturach cienkościennych. Metoda ta oparta jest na założeniu 
jednowymiarowego przewodzenia ciepła w elementach struktur. Omówiono 
sposoby rozwiązania zagadnienia dla dwóch wariantów: w wariancie A pojem- 
ność cieplna czynników omywających elementy struktur była nieskończenie 
wielka, w wariancie B pojemność czynników mogła być ograniczona. Dla 
wybranej struktury w wariancie A dokonano porównania rozwiązania uprosz- 
czonego jednowymiarowego z rozwiązaniem dokładnym dwuwymiarowym; na 
tej podstawie podano zakres stosowania rozwiązań jednowymiarowych. 
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W przypadku gdy liczba Biota liczona względem czynnika grzejnego 
Bi, < 0,5 i gdy stosunek liczb Biota względem czynnika grzejnego i grzanego 
Bi,/Bi > 20 bądź Bi < 0,01, maksymalne różnice między wartościami tem- 
peratury określonymi na powierzchni górnej i dolnej oraz wartością określoną 
z modelu jednowymiarowego względem temperatury odniesienia mogą wyno- 
sić około 2,0%. 

Należy podkreślić, że stosunkowo prosto można wyznaczyć pole tem- 
peratury dla struktury złożonej z małej liczby elementów, natomiast dla 
struktur o większej ich liczbie zagadnienie się komplikuje. Z tego powodu 
w strukturze wyodrębnia się elementy powtarzalne i dla nich wykonuje 
obliczenia. | 


8. WYKAZ OZNACZEŃ 


— wymiar liniowy, współczynnik 

— wymiar liniowy zredukowany, współczynnik 
— macierz współczynników równania 

— długość pasma, współczynnik 

— wymiar liniowy zredukowany 

— szerokość elementu, współczynnik 

— wymiar liniowy zredukowany 

— grubość elementu, współczynnik 

— współczynnik 

— wektor kolumnowy wyrazów wolnych 
— stała 

— współczynnik 

— liczba naturalna 

— |liczba naturalna 

— liczba naturalna, funkcja 

— wymiar liniowy 

liczba naturalna 

— |liczba naturalna 

— funkcja 

— strumień cieplny 

— strumień cieplny 

— temperatura czynnika 

— temperatura płyty 

— wektor kolumnowy niewiadomych temperatur 
— temperatura Średnia 

— pojemność cieplna strumienia czynnika 
— współrzędna wzdłuż szerokości płyty 
— współrzędna zredukowana 

— współrzędna wzdłuż długości płyty 

— współrzędna zredukowana 


<= HRZJJJTOSEO SZR" "O MZUROASTEWSNANS 
| 
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— współrzędna wzdłuż grubości płyty 
— współrzędna zredukowana 

— współczynnik wnikania ciepła 

— współczynnik 

współczynnik 

— współczynnik przewodności cieplnej 
— współczynnik 

— temperatura zredukowana 


QE RMUMWLXLNN 
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8.1. MODUŁY BEZWYMIAROWE 


| 
Bi= > — liczba Biota. 


A 
Indeksy dotyczą: g — grzania, i — numeru węzła i — kolejnego 
elementu (czynnika), j — numeru węzła, ij — pasma zawartego pomiędzy 
węzłami ij, k — liczby czynników graniczących z elementem, m — liczby 
czynników wpływających na temperaturę pasma, o — otoczenia, | — pierw- 
szego elementu (czynnika), 2 — drugiego elementu (czynnika). 


WŁADYSŁAW NOWAK, WŁADYSŁAW SZAFLIK 


A SIMPLIFIED METHOD OF SPECIFYING THERMAL FIELDS 
IN THIN-WALLED BAND STRUCTURES 


Summary 


While heating containers' panels, bottoms and walls of various kinds a method of heating 
medium guidance in parallel semi-cylindrical or angular profiles, welded to one of the surfaces of 
the heated elements, is employed. 

In the solutions listed above, specifying exact amounts of transmitted heat is generally very 
complicated, and in some cases virtually impossible. Therefore, in complex structures, made from 
reproducible thin-walled elements, it is possible to use approximate methods taking only 
unidirectional thermal conduction into consideration. Band structures are such thin-walled 
structures in which band-shaped elements can be distinguished differing in heat exchange 
conditions. 

A simplified analytical method of thermal field determination in thin-walled band structures is 
presented in the paper. The method is based on the assumption of unidirectional thermal 
conduction in the elements of the structure. The methods of the problem solving for two variants 
were discussed: in variant A — thermal capacity of agents washing the structures elements was 
infinitely high; in variant B thermal capacity of agent flows washing the structures element had 
a finite value. In variant A the comparison between the simplified unidirectional solution and 
non-simplified bi-directional one for a chosen structure was made. 

The analysis allowed to determine the scope of application of unidirectional solutions, 
for which maximum differences between the mean value, calculated with assumption of 
a bi-directional thermal field, and the value determined by the unidirectional model can be of the 
order of 2.0%. 
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BUAĄBICJIAB HOBAK, BJIAĄBCJIAB IIAQJIAK 


COKPAIIEHHBIH METOĄ OIIPEJEJIEHMA IIOJLA TEMIIEPATYP 
B IIOJIOCATBIX TOHKOCTEHHBIX CTPYKTYPAX 


Pe3oMe 


IIpu nonorpeBaHuH pa3HOTO BHJIA IIJIHT, JĄHHIII H CTEHOK pe3epByapoB IIpHMEHAETCA CIIOCOÓ 
BBEJIEHHA TEIIJIOHOCHTEJIA B NapaJLIEJIbBHBIX IOJIYNMJIAH1DHYECKHX HJIA yTJIOBbIX, HPpHBADEHHBIX 
K OJIHOŃ H3 IIOBEpxHOCTEeŃ IONOTPEBAEMBIX 3JIEMEHTOB, HpoQH1Ax. B Bblliieyka3aHHBIX PEILIEHHAX 
TOUHO€ OIIpEIEJIEHHE KOJIMUECTBA IIEpeNaBaeMOTO TEILIA ABJIAETCA B OÓLIEM CJIy4ae OUEHB 
CJIOXKHBIM, A HHOTJIA BOOÓLINE HEBO3MOXHBIM. IIO3TOMy HpH CJIOXKHbIX KOHCTDYKLUH4AX, CHEJIAHHBIX 
H3 IIOBTODHBIX TOHKOCTEHHbIX 3JIEMEHTOB, MOXHO IIDHMEHHTb HpHÓJIMXeHHble METOJIBI, YHHTbIBa- 
IOLIKE TENJIO TOJIbBKO B OHHOM HalpaBJIeHHH. [IOJOCATbie CTPYKTYPKI — 3TO TAKHE CTPYKTYDBI, 
B KOTODPBIX MOXHO BbI/EJIMTb J3JIEMEHTBlI B BHJIE IIOJIOC, OTJIAUAŁOLNKECA YCJIOBHAMH TEIIJIOOOMEHA. 

B paóoTe npeĄcTaBJIeH COKPALIeHHbBIA AHAJJATHUECKHŃ METOJĄ ONDPEHNEJJEHAA IOJIA TEM- 
repaTyp B IIOJIOCATBIX TOHKOCTEHHbIX CTDYKTYpaX. DTOT METOJĄ OCHOBAH Ha IHpENIIOJOXEHHK, UTO 
TEILIO B 3JIEMEHTAX CTPYKTYDPBI HPOXOJIKT B OJHOM HalpaBJIeHHH. OTOBODpEHBI CIIOCOÓBI peLIIeHHA 
HpOÓJIEMbI NIA HAByX BApHAHTOB: B BADHaHTe A — TEeIIOBOŃ OÓbeM (haKTOPOB, OMBIBAIOIIIHX 
3JIEMEHTbI CTPYKTyp, ÓbiI Oe3TDAaHHUHO ÓOJIBIIOŃ, B BapuaHTe B TENIOBOŃ OÓbeM IMOTOKA 
(baKTOPOB, OMBIBAFOLIHX 3JIEMEHTBI CTDYKTyp, HMeeT OTDPAHHUEHHOE 3HAUeHKHE. JIJIA BapHaHTa A, 
nia u30paHHoń CTPykTypbl IpENIOXKEHO CpAaBH€HHE COKPAIIIEHHOTO OJ]HOMEpHOTO peLNieHHA 
C IIOJIHBIM J1ByMEpHBIM DpeLIIEHHEM. 

Ha OCHOBe aHaJIH3a ONpENEJIEHA OÓJIACTb IDUMEHEHHA ONHOMEpHOTO peLIeHHA, NJIA KOTOPOŃ 
MAaKCHMAJIbHAA pAa3HHIIA MEXJIy CpPEHHHM BBIUMCJIEHHbBIM 3HA4CHHEM IIDH IHpENIOJIOXKEHHH JĄBY- 
MepHOTO IIOJI4 TEMIIepaTyp M OIIpEIEJIEHHOTO 3HAa4CHHA H3 ONHOMEpHOŃ MOJIEJIM, MOXKET ObITb 
nopaaąka 2,0%. 
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Streszczenie. Przedstawiono zagadnienie dekompozycji problemu poliop- 
tymalizacji dyskretnej. Dekompozycji mogą podlegać następujące elementy procesu 
optymalizacji: wektor zmiennych decyzyjnych, wektor funkcji celu, obszar rozwiązań 
dopuszczalnych, metoda optymalizacji oraz obiekt optymalizacji. Zastosowanie 
w jednym zadaniu różnych form dekompozycji znacznie przyśpiesza uzyskanie 
rozwiązania. W pracy zestawiono 112 kryteriów możliwych do zastosowania 
w optymalizacji konstrukcji i ich elementów. Podzielono je na sześć grup: nadrzęd- 
ne — 9, konstrukcyjne — 35, technologiczne — 16, ekonomiczne — 22, ekolo- 
giczne — 13 i funkcjonalne — 20. Pokazano możliwości dekompozycji zadania na 
przykładzie hali biurowo-magazynowej. 


Słowa kluczowe: polioptymalizacja — dekompozycja — kryterium — dys- 
kretny — niezdominowany — konstrukcja. 


1. WPROWADZENIE 


Tradycyjne projektowanie konstrukcji budowlanych ugruntowało stosowa- 
nie dekompozycji obiektu na współpracujące ze sobą elementy, np. szkielet — 
stropy — ściany — fundamenty, i analizy każdego z nich osobno. Projek- 
towanie to uporządkowane zostało sposobem przekazywania obciążeń: z naj- 
wyższych kondygnacji, poprzez niższe na fundamenty lub z przęseł na podpory. 
Dzięki podziałowi obiektu na elementy możliwe stało się równoległe projek- 
towanie [7] kilku elementów jednocześnie. Proces taki przy odpowiedniej 
koordynacji poszczególnych zadań lokalnych znacznie upraszcza analizę oraz 
przyspiesza opracowanie końcowego projektu globalnego. 

W projektowaniu optymalnym [4] dąży się do uzyskania takiej konstruk- 
cji. która oprócz spełnienia podstawowych przepisów normatywnych jest 
najlepsza ze względu na ustalone kryterium, np. minimum kosztu. W projek- 
towaniu polioptymalnym [3, 28] poszukujemy rozwiązań niezdominowanych, 
czyli możliwie najlepszych ze względu na kilka konfliktowych między sobą 
kryteriów, np. minimum masy i maksimum sztywności. Spośród nich za 
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pomocą dodatkowego kryterium wybierane jest rozwiązanie preferowane, 
którego realizacja jest najbardziej racjonalna. 

Coraz częściej spotyka się postulat, aby zadanie optymalizacji konstrukcji 
budowlanych formułowane było jako dyskretne i wielokryterialne [ 7, 11, 17, 24, 
31]. W przypadku analizy obiektu liczba istotnych zmiennych decyzyjnych 
powoduje, że w obszarze dopuszczalnym znajduje się od kilku tysięcy do kilku 
milionów rozwiązań. Przegląd i ocena wszystkich jest najczęściej niemożliwa 
przy dzisiejszym stanie techniki obliczeniowej. Jednym ze sposobów rozwiązania 
problemu jest, podobnie jak w tradycyjnym projektowaniu, dekompozycja 
zadania, analiza kilku mniejszych zadań lokalnych i w końcu ich koordynacja 
w taki sposób, aby uzyskane rozwiązania były poprawnym wynikiem zadania 
globalnego. Powstałe w ostatnich latach metody dekompozycji dotyczą głównie 
optymalizacji jednokryterialnej [2, 4, 5, 9, 10, 13, 14, 16, 17, 32—34]. 
W zagadnieniach wielokryterialnych [5, 6, 8, 12, 15, 19, 21, 23, 26, 28—31] 
problem dekompozycji jest bardziej złóżony, ponieważ obiektywnym wynikiem 
rozwiązania zadań lokalnych jest zbiór rozwiązań niezdominowanych, który 
należałoby uwzględniać przy optymalizacji globalnej [6, 23, 26]. Ponadto 
zmienne decyzyjne i funkcje celu zadań lokalnych nie muszą się pokrywać ze 
zmiennymi decyzyjnymi i funkcjami celu zadania globalnego [28]. Dekom- 
pozycji podlegać mogą wszystkie elementy procesu optymalizacji: wektor 
zmiennych decyzyjnych [8, 11, 12, 16, 29], wektor funkcji celu [8, 12, 29, 34], 
obszar dopuszczalny [8, 14, 26], metoda optymalizacji [3, 7, 10] oraz obiekt 
optymalizacji [2, 8, 9, 14, 15, 17, 23]. W praktyce wykorzystuje się często kilka 
form dekompozycji w jednym zadaniu, nie rozgraniczając ich [7, 10]. Prowadzi 
to niekiedy do opracowania oryginalnych algorytmów, przeznaczonych do 
rozwiązania wąskiej klasy zagadnień. Szerszy przegląd koncepcji dekompozycji 
stosowanych w problemach optymalizacji zawierają prace [6, 21, 26, 29]. 


2. SFORMUŁOWANIE ZADANIA POLIOPTYMALIZACJI 


Na podejmowane decyzje wpływ ma zazwyczaj wiele czynników. Decyden- 
towi zależy na tym, aby zaakceptowany przez niego projekt przyniósł jak 
największe efekty, np. ekonomiczne, bądź też aby jak najbardziej łagodził 
niekorzystne, np. ekologiczne, skutki podjęcia określonych inwestycji. Przygoto- 
wanie odpowiednich danych, ułatwiających podjęcie trafnej decyzji, może być 
przeprowadzone w wyniku rozwiązania odpowiednio zdefiniowanego zadania 
optymalizacji wielokryterialnej. Zadanie takie jest sformułowane, jeżeli dla 
danego obiektu (np. budynku, mostu, hali) określone zostaną zmienne decyzyj- 
ne, ograniczenia 1 kryteria optymalizacji oraz przyjęta będzie relacja porząd- 
kująca oceny rozwiązań [1]. 

Zmienne decyzyjne są wielkościami definiującymi analizowany obiekt, 
podlegającymi kształtowaniu w procesie optymalizacji. Zmienne decyzyjne mają 
postać wektora 


x = [X1 X2, DEETE <a Fędl la 
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w N-wymiarowej przestrzeni zmiennych decyzyjnych A c R”. Każdemu punk- 
towi x w przestrzeni rozwiązań A odpowiada obiekt opisany wektorem 
N zmiennych decyzyjnych (ryc. I, N=2). Obszar dopuszczalny X stanowi 
część przestrzeni A (X c 4), wyznaczoną przez ograniczenia nałożone na 
zmienne decyzyjne (ryc. 1). Ograniczenia mają postać nierówności 


gwixitś 0; gdzie" JEŻ, =NM 2a l zas 
i równości 

Gł) =D. gdzie. W EŚ, = 3, lono TJ 
Funkcja celu f(x) jest funkcją opisującą pewną właściwość obiektu op- 
tymalizacji, przyporządkowującą poszczególnym stopniom realizacji danego 
celu wartości stanowiące podstawę oceny rozpatrywanego obiektu. Optymali- 
zacja według określonego kryterium to poszukiwanie minimum lub maksimum 
funkcji celu f(x) w zbiorze decyzji dopuszczalnych X. 


X2 f(x) B 
y =f(x) o y2€ ŻA 
yrtA 
Y1€ > 
m(x)=0 YcB 
V; 
min f(x) o 
A 
X1 min f,(x) f1(X) 


Ryc. 1. Graficzna interpretacja zadania optymalizacji dwukryterialnej (J = 2) z dwoma zmiennymi 
decyzyjnymi (N = 2) 


Jeżeli wymagana jest ocena obiektu ze względu na J > 1 kryteriów, to 
zadanie nosi nazwę optymalizacji wielokryterialnej (polioptymalizacji) i jest 
traktowane jak zadanie optymalizacji J-wektorowej funkcji celu 


10 =L10), .... 50%), ....5,(9]'. 


J-wymiarowa przestrzeń, w której odwzorowane są wartości funkcji celu f(x), 
jest przestrzenią celu B = R”. Zbiór ocen decyzji Y jest tą częścią przestrzeni 
celu B, w której położone są wartości funkcji celu odpowiadające dopuszczal- 
nym wartościom zmiennych decyzyjnych (ryc. 1, J = 2). Odwzorowanie zbioru 
X na zbiór Y jest suriekcją, co znaczy, że przeciwobraz każdego punktu 
y składa się co najmniej z jednego punktu należącego do zbioru X: 


fiX—>Y, x>y=f(x, xeXcACGR" yeYcBcR. 


Oceny decyzji y zawarte w zbiorze ocen Y podlegają porównaniu między sobą 
za pomocą relacji % porządkującej zbiór Y [1]. Relacja ta powinna umożliwić 
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porównanie dwóch dowolnych ocen i orzeczenie, która z nich jest dominująca, 
a w przypadku ich indyferentności stwierdzenie, że oceny są niezdominowane. 
Najczęściej przyjmowana jest stożkowa relacja porządkująca %, = <, ze 
stożkiem A określonym przez orthant przestrzeni R” [1] 

AZ MAN2, SAPER 20: aIbO: ARAB NEGO 


J 
Zbiór ocen niezdominowanych (optymalnych w sensie Pareto) Yyp, stanowiący 
obiektywny wynik zadania polioptymalizacji, jest wtedy położony na brzegu 
zbioru ocen Y (rys. 1): 
Yb =J/(Xnp) = Opt/(x) = 


xEX 
= (yNpEY: Idy, EYAy; 4 yNp, takie Że yyp€Ey, +4]. 


Zbiorowi Yyp odpowiada w przestrzeni decyzji A zbiór decyzji niezdominowa- 
nych Xpp (ryc. 1): 


Xin =f" (Wyp) = (XNpPEX: YNb = J(XNp)€ YNbj. 


Zbiór Pareto Yyp zawiera najczęściej wiele ocen decyzji. Spośród nich na 
podstawie dodatkowego kryterium, np. minimalnej odległości od oceny ideal- 


nej Vid> 
Mid rej LV1.1; V2.2> ..., VIu], 
Via,j = Vij = OpttJ,(x): xeX, je f), Opte(Max, Min), 


wybierana jest ocena preferowana y, i odpowiadająca jej w zbiorze X decyzja 
preferowana x. 


3. DEKOMPOZYCJA WEKTORA ZMIENNYCH DECYZYJNYCH 


W dekompozycji wektora zmiennych decyzyjnych wykorzystuje się właś- 
ciwość monotoniczności zmiennych [29]. Załóżmy (ryc. 2), że mamy zadanie 
polioptymalizacji z N> 1 zmiennymi decyzyjnymi. Przez Yyp(n) oznaczymy 
zbiór ocen niezdominowanych, a przez X yp(n) zbiór rozwiązań niezdominowa- 
nych dla n < N zmiennych decyzyjnych (ryc. 3). Przestrzeń zmiennych decyzyj- 
nych A ograniczymy w ten sposób, że spośród N składowych wektora 
x ustalamy dowolną, np. (N—1), składową. Ograniczenie to powoduje, że 
N-wymiarowy obszar dopuszczalny X(N) jest rzutowany na przestrzeń 
(N—1)-wymiarową, czyli powstaje nowy obszar dopuszczalny X (N—1). Roz- 
wiązujemy teraz zadanie polioptymalizacji dla (N— 1) zmiennych decyzyjnych, 
tzn. określamy zbiory Typ(N— 1) i Xyp(N—1). Z uwagi na warunek porządku 
liniowego wprowadzony przez stożkową relację < „ z ograniczonego obszaru 
dopuszczalnego X (N— 1) w żaden sposób nie można poprawić całego zbioru 
Yyp(N). Najczęściej sytuacja pogarsza się lub dla niektórych ocen niezdomino- 
wanych pozostaje bez zmian. 
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Ryc. 2. Dekompozycja wektora zmiennych decyzyjnych w zadaniu jednokryterialnym 


Rozwiązując zadanie optymalizacji postępujemy odwrotnie, zaczynając od 
jednej lub dwóch zmiennych decyzyjnych, a następnie dodając kolejne. 
Dołączenie kolejnej (n + 1) zmiennej decyzyjnej, gdzie (n-+- 1) < N, może jedynie 
poprawić (tzn. zmniejszyć przy minimalizacji, a zwiększyć przy maksymalizacji) 
wartości funkcji celu. 

W przypadku zadań wielokryterialnych dołączenie kolejnej zmiennej decy- 
zyjnej poprawia zbiór ocen niezdominowanych (ryc. 3a, c). Jeżeli nowa zmienna 
decyzyjna jest nieaktywna, tzn. nie generuje lepszych wartości funkcji celu, to 
zbiór Yyp(n) nie zmienia się (ryc. 3b). 


fz(x) f(x) 


Y(n+1) 


f,00) f,(X) 


Ryc. 3. Zbiory ocen rozwiązań dla n i (n+1) zmiennych decyzyjnych 


4. DEKOMPOZYCJA WEKTORA FUNKCJI CELU 


Właściwość monotoniczności wektora funkcji celu polega na przynależ- 
ności do zbioru ocen niezdominowanych Yyp poszczególnych niezdomino- 
wanych zbiorów cząstkowych Yyp(j), j <J, gdzie J jest liczbą kryteriów 
zadania [1, 29]. Symbolicznie właściwość monotoniczności celów można 
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zapisać następująco: 
VJ<J, Pol) = do(i+1) © Typ. 


Zależność pomiędzy zbiorami ocen niezdominowanych o różnej liczbie kryte- 
riów polega na tym, że przy rzutowaniu zbioru Yyp(f+1) na przestrzeń R” 
zbiór Yyp(j) jest zbiorem ocen niezdominowanych rzutowanego zbioru o więk- 
szej wymiarowości. Na rycinie 4 pokazano zbiór ocen niezdominowanych 
w przestrzeni R* oraz niezdominowane zbiory ocen cząstkowych Yyp(2), 
Yxp (1). Indeksami górnymi oznaczono numery kryteriów cząstkowych: 


Typ (2) e fYN5 "> VND» YND *). 


fz(x) 


P Ś? (x) 


Ryc. 4. Pełne i cząstkowe zbiory ocen niezdominowanych zadania trójkryterialnego 


W przypadku granicznym j = 1, elementy zbioru Yyp(1) stanowią indywidual- 
ne wartości ekstremalizujące wybrane funkcje celu: 


Yyp (1) = (YND; YAD; YWp) = RZE y2, y3). 


Praktyczne wykorzystanie dekompozycji wektora funkcji celu polega na 
tym, że w pierwszej kolejności określa się ocenę narożną Yyp = yi, a następnie 
cząstkowy zbiór niezdominowany Yjp *. Później wyznacza się kolejne cząst- 
kowe zbiory niezdominowane dwukryterialne: Yyp*, Ysp *, Ib", VNb*, 
YJp * itd. (ryc. 5a). Zazwyczaj określenie tych zbiorów satysfakcjonuje prowa- 
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dzącego analizę optymalizacyjną i rezygnuje on ze znalezienia zbiorów 
niezdominowanych trój- i więcej kryterialnych (np. Yyp *"*, ryc. 5b), będących 
podzbiorami ostatecznego zbioru ocen niezdominowanych Yyp = Typ(J). 


Ryc. 5. Cząstkowe zbiory ocen niezdominowanych zadania czterokryterialnego (a) i trójkryterial- 
nego (b) 


W tabeli 1 zestawiono kryteria możliwe do zastosowania w optymalizacji 
konstrukcji budowlanych [27]. Podzielono je na sześć grup: nadrzędne [34], 
konstrukcyjne, technologiczne, ekonomiczne, ekologiczne i funkcjonalne. Kry- 
teria nadrzędne określone są we wstępnym etapie formułowania zadania 
optymalizacji i wyznaczają pewne ogólne cele, jakie powinien spełniać rozpat- 
rywany obiekt. Kryteria te przyjmują na ogół formę werbalną, a po sprecyzo- 
waniu stają się zadaniowymi kryteriami optymalizacji. Przykładowo minimum 
oczekiwanego kosztu obiektu, po szczegółowej analizie rynku, może przyjąć 
formę minimum kosztu realizacji obiektu i minimum zapotrzebowania na 
energię przy eksploatacji. W przypadku rozwiązywania problemów lokalnych, 
takich jak określenie typu mostu czy schematu statycznego budynku, potrzeb- 
ne są czasem indywidualne kryteria nadrzędne, które mogą być wybrane 
z kryteriów zadaniowych pokazanych w tabeli 1. 

Od obiektu optymalizacjj wymaga się najczęściej, aby spełniał wiele 
sprzecznych ze sobą kryteriów. Zadaniowe kryteria optymalizacji, definiowane 
w wektorowej funkcji celu f(x), pozwalają na formułowanie problemów 
wielokryterialnych. W przypadku zadań z zakresu budownictwa wyróżniono 
pięć grup kryteriów (tab. 1). Ich znaczna liczba powoduje, że nie mogą być one 
wszystkie jednocześnie brane pod uwagę. Odpowiedni dobór od 3 do 7 kryte- 
riów, najlepiej z różnych grup, gwarantuje wszechstronną ocenę obiektu [27]. 
Stosowanie większej liczby kryteriów oceny powoduje nadmierne zwiększenie 
liczebności rozwiązań niezdominowanych, przez co wybór rozwiązania prefe- 
rowanego jest niekiedy znacznie utrudniony [34]. Istotne z punktu widzenia 


Tabela 1. Kryteria optymalizacji konstrukcji budowlanych 


Kryteria nadrzędne 


U 


| Zadaniowe kryteria optymalizacji ] 
4 
Ekonomiczne i Funkcjonalne M Konstrukcyjne | | Technologiczne | | Ekologiczne 


Kryteria nadrzędne 


1) max spełnienia potrzeb ludności i gospodarki w zakresie infrastruktury technicznej 
i społecznej 

2) min oczekiwanego kosztu obiektu (planowanie, zakup gruntu, projektowanie, realizacja, 
eksploatacja, zmiany funkcji, ew. awaria lub katastrofa, rozbiórka, odzysk materiałowy, 
utylizacja) 

3) max preferencji dla technologii o zamkniętym obiegu oraz ograniczających emisję szkod- 
liwych substancji stałych, pyłów i gazów 


4) max bezpieczeństwa lub niezawodności konstrukcji przy ustalonych nakładach 

5) min wskaźnika efektywności energetycznej obiektu lub min energochłonności obiektu 
6) max przystosowania obiektu do potrzeb człowieka i procesów technologicznych 

7) min ingerencji w lokalny i globalny ekosystem w czasie realizacji i eksploatacji obiektu 
) 


8) max zysków (np. finansowych) i min strat (np. ekologicznych) wynikających z lokalizacji | 


obiektu 
9) max wzrostu ekonomicznego i konkurencyjności krajowych firm budowlanych o zróż- 
nicowanej wielkości 


Kryteria ekonomiczne Kryteria funkcjonalne 


1) max użyteczności przy ustalonych nakła- 1) max przystosowania do psychofizycznych 
dach i estetycznych potrzeb człowieka 

2) min nakładów przy ustalonych wymaga- 2) max dopasowania do pełnionych funkcji 
niach użyteczności 3) max łatwości eksploatacji 

3) max efektywności ekonomicznej inwesty- 4) max komfortu cieplno-wilgotnościowego 
cji i akustycznego pomieszczeń 

4) min kosztów dokumentacji projektowej 5) max bezpieczeństwa ruchu pieszego 

5) min kosztów materiału (przeciwpoślizgowość, bezkolizyjność) 

6) min kosztów zaopatrzenia materiałowego 6) max wytłumienia hałasu 

7) min kosztów transportu i montażu ele- 7) max wytłumienia oddziaływań dynamicz- 
mentów nych 

8) min oczekiwanych kosztów realizacji 8) max oświetlenia naturalnego wewnątrz 
obiektu obiektu 

9) min oczekiwanych kosztów eksploatacji 9) max dostosowania trwałości elementów 
obiektu do czasu eksploatacji obiektu 


| 10) min kosztów energii potrzebnej do ogrza- | 10) max przystosowania do potrzeb osób 


nia, wentylacji i oświetlenia obiektu niepełnosprawnych 
11) min kosztów ewentualnej awarii lub kata- | 11) max elastyczności funkcji obiektu 
strofy 12) max bezpieczeństwa przeciwpożarowego 
12) min nakładów na konserwacje i remonty | 13) max estetyki i wrażenia architektonicz- 


| 13) min kosztu utylizacji obiektu nego 


| 
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Tabela 1 (cd.) 


14) min kosztu zmian funkcjonalnych w obie- 
kcie eksploatowanym 

15) max wykorzystania potencjału firmy bu- 
dowlanej 

16) min czasu realizacji obiektu 

17) max sprawności rozwiązań organizacji 
inwestycji 

18) max wydajności pracy przy realizacji 
obiektu 

19) max możliwości kompleksowego wyko- 
nawstwa 

20) max zysku firmy budowlanej 

21) max odzysku materiałów i elementów 
z demontażu konstrukcji 

22) min kosztów dojazdu z osiedla miesz- 
kaniowego do dzielnic przemysłowych 
i biurowych 


14) max spełnienia wymogów norm jakości 

15) max trwałości barw 

16) max odporności na zabrudzenia i zapy- 
lenia 

17) max powierzchni użytkowej netto 

18) max przepustowości arterii komunikacyj- 
nych 

19) min odległości osiedli mieszkaniowych od 
zakładów pracy 

20) min czasu i max komfortu dojazdu do 


pracy 


Kryteria konstrukcyjne Kryteria technologiczne 


1) min masy lub objętości materiału 

2) max stateczności konstrukcji i elementów 

3) max obciążenia krytycznego przy ustalo- 
nej objętości materiału 

4) max wyrównania stateczności elementów 

5) min wrażliwości na zmiany parametrów 
konstrukcyjno-materiałowych 

6) max sztywności, min podatności lub max 
momentów bezwładności przekrojów 

7) max wyrównania naprężeń lub wytężeń 

8) max wyrównania naprężeń brzegowych 

9) max wyrównania odkształceń jednostko- 
wych 

10) max wyrównania potencjałów jednostko- 
wych 

11) max obniżenia koncentracji naprężeń 

12) min energii sprężystej (min potencjału) 

13) min przemieszczeń wybranych punktów 
konstrukcji 

14) max niezawodności wzmocnienia pod- 
łoża 

15) max nośności sprężysto-plastycznej 

16) max przystosowania do obciążeń zmien- 
nych 

17) max pierwszej częstości drgań własnych 

18) min objętości lub masy przy określonej 
częstości podstawowej 

19) max różnicy między pierwszą częstością 
drgań własnych i częstością obciążenia 
wymuszającego 


1) max technologiczności konstrukcji (po- 
wtarzalność, modułowość, tolerancje wy- 
miarowe, łatwość transportu i montażu, 
dostosowanie do katalogów wyrobów) 

2) max trwałości obiektu 

3) max typizacji elementów i obiektów 

4) min zróżnicowania przekrojów elemen- 
tów 

5) min liczby połączeń lub liczby elementów 

6) min liczby połączeń montażowych 

7) max wykorzystania surowców lub mini- 
malizacja odpadów materiałowych 

8) max dostosowania produkcji budowlanej 
do międzynarodowych norm jakości ISO 

9) max odporności na niefachowe wyko- 
nawstwo i eksploatację 

10) max odporności na niekorzystne wpływy 
klimatyczne (wodoodporność, mrozood- 
porność, odporność na korozję, wpływy 
termiczne, działanie wilgoci) 

11) max odporności na agresywne środowisko 

12) max odporności na działanie podwyższo- 
nej lub obniżonej temperatury 

13) max odporności na Ścieranie 

14) max odporności na starzenie się i zużycie 

15) max odporności na zarysowania 

16) max zdolności do wykonywania napraw 
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Tabela 1 (cd.) 


| Kryteria konstrukcyjne Kryteria ekologiczne 


20) max różnicy między zadaną częstością 
siły wymuszającej a sąsiednimi częstoś- 
ciami drgań własnych 

21) min amplitudy drgań lub max wytłumie- 
nia drgań 

22) max odporności na nierównomierne osia- 
dania 

23) max dostosowania do obciążeń sejsmicz- 
nych i parasejsmicznych 

| 24) max odporności na odkształcenia plas- 
tyczne 

25) max odporności na zarysowania 

26) min rozwarcia rys elementów żelbeto- 
wych 

27) min mimośrodów w przekazywaniu ob- 
ciążeń 

28) max dopasowania typu i kształtu kon- 
strukcji do wielu stanów przenoszonych 
obciążeń 

29) min wrażliwości na oddziaływania ter- 
miczne 

30) min wrażliwości na oddziaływania reolo- 
giczne 

31) max odporności ogniowej 

32) max odporności na działania mechanicz- 
ne 

33) max harmonizacji przepisów projektowa- 
nia 

34) max efektywności wzmocnienia konstru- 
kcji lub jej elementu 

35) min funkcji uchybu identyfikacji obiektu 
(np. min sumy odchyleń kwadratowych 
pomiędzy odkształceniami modelu a od- 
kształceniami konstrukcji od obciążenia 
próbnego) 


1) min zapotrzebowania na energię przy 
eksploatacji obiektu (ogrzewanie lub 
chłodzenie) 

2) min energochłonności realizacji obiektu 

3) max zużycia surowców wtórnych, odpa- 
dowych, odnawialnych, lokalnych 

4) min stosowania materiałów szkodliwych 
dla zdrowia człowieka oraz dla Środo- 
wiska 

5) max wykorzystania procesów technologi- 
cznych przyjaznych dla środowiska natu- 
ralnego 

6) min odpadów materiałowych zaśmiecają- 
cych środowisko 1 nie ulegających bio- 
degradacji 

7) min powierzchni placu budowy 

8) max przystosowania do istniejącego oto- 
czenia naturalnego 

9) max wykorzystania ekologicznych środ- 
ków transportu 

10) max sprawności urządzeń oczyszczają- 
cych 

11) max szerokości stref ochronnych 

12) min ingerencji w strefy chronione 

13) min wpływu na środowisko obiektów 
tymczasowych towarzyszących inwestycji 


danego obiektu kryteria oceny, lecz takie, które nie weszły w skład zadaniowe- 
go wektora funkcji celu f(x), przechodzą do ograniczeń zadania. Przykładowo 
często rezygnuje się z kryterium maksimum sztywności, włączając do ograni- 
czeń warunek nieprzekroczenia przemieszczenia granicznego. 


5. DEKOMPOZYCJA OBSZARU DOPUSZCZALNEGO 


Dekompozycja obszaru dopuszczalnego wykorzystywana jest w sposób 
naturalny w przypadkach, gdy jest on wielospójny, to znaczy złożony z kilku 
rozłącznych części. Poszukiwania prowadzi się wtedy niezależnie w poszczegól- 
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nych podobszarach obszaru dopuszczalnego, a w końcowej fazie wybiera się 
rozwiązanie globalne x, (ryc. 6a). 

Jeżeli rozwiązujemy zadanie polioptymalizacji, to w wyniku otrzymujemy 
podzbiory ocen odpowiadające poszczególnym podobszarom dopuszczalnym. 
Ostatecznie wybierany jest taki zbiór rozwiązań niezdominowanych, którego 
oceny tworzą brzeg Pareto zbioru Y (ryc. 6c). Mogą to być oceny od- 
powiadające rozwiązaniom położonym w jednym lub kilku podobszarach 
dopuszczalnych. 


X2 Xx | 


= 


(00/0822 


f(x) X1 | Xi 09 


Ryc. 6. Dekompozycja obszaru dopuszczalnego wielospójnego dla zadania jednokryterialnego (a) 
i wielokryterialnego (b, c) 


Dekompozycję obszaru dopuszczalnego należy przeprowadzić w przypad- 
kach, gdy mamy informacje mogące świadczyć o przewężeniach (ryc. 7). Są one 
efektem zastosowania pewnych ograniczeń, np. technologicznych. Przewężenia 


x2 


Ryc. 7. Przewężenie obszaru dopuszczalnego w zadaniu jednopoziomowym (a) i wielopoziomo- 
wym (b) 
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powstają również w zadaniach wielopoziomowych, w których poszczególne 
elementy konstrukcji opisane są różnymi zmiennymi decyzyjnymi lokalnymi. 
Rozczłonkowany obszar dopuszczalny połączony jest dzięki zmiennym decy- 
zyjnym globalnym. Dekompozycję przeprowadza się w przewężeniach między 
zadaniami lokalnymi a zadaniem globalnym. 

Inną koncepćję dekompozycji obszaru dopuszczalnego można zastosować, 
wykorzystując niektóre standardowe metody optymalizacji, np. modyfikowaną 
metodę Monte Carlo. W obszar dopuszczalny wpisana zostaje regularna siatka 
dzieląca cały zbiór X na M N-wymiarowych kostek (ryc. 8). W każdej kostce 
K,„,me.mM =1,2,...,m,..., M), przeprowadza się losowanie. Do następ- 
nego etapu optymalizacji wybrana zostaje kostka, w której funkcja celu 
osiągnęła najmniejszą wartość (ryc. 8a). W przypadku polioptymalizacji 
wybierane zostają te kostki, w której rozwiązania uzyskały oceny niezdomino- 
wane (ryc. 8b). W drugim etapie losowanie przeprowadza się w obszarze 
ograniczonym do jednej lub kilku kostek, w zależności od liczby kryteriów 
optymalizacji. Kolejne etapy polegają na zmniejszaniu wymiarów kostek 
i przesuwaniu ich w kierunku rozwiązania zadania. W przypadku problemu 
dyskretnego warunkiem stopu jest osiągnięcie wymiarów kostki odpowiadają- 
cych otoczeniu danego punktu. 


Ryc. 8. Podział obszaru dopuszczalnego na MN-wymiarowych kostek w zadaniu jednokryterial- 
nym (a) i wielokryterialnym (b) 


6. DEKOMPOZYCJA METODY OPTYMALIZACJI 


Dekompozycja metody optymalizacji polega na współdziałaniu kilku 
strategii przeszukiwania obszaru dopuszczalnego połączonych w jednym 
algorytmie optymalizacyjnym [11, 25, 30, 36, 37]. Przykładowo zastosowanie 
metody Monte Carlo do wyznaczenia dobrego punktu startowego, meto- 
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dy Gaussa-Seidela do wyznaczenia rozwiązania znajdującego się w pobliżu 
optimum oraz metody siatki przeszukującej otoczenie danego punktu 
w celu uściślenia wyniku może stanowić efektywne narzędzie optymaliza- 
cyjne [11]. 

Tworzone w ostatnich latach metody optymalizacji, jak np. metoda 
ortogonalno-diagonalna [25 |, mają szeroko rozbudowane sterowanie, umoż- 
liwiające uruchamianie mniej lub bardziej zaawansowanych sposobów prze- 
szukiwania obszaru dopuszczalnego, w celu wyznaczenia kierunku 'poprawy 
funkcji celu. Sterowaniem tym nadzoruje system ekspercki, wspomagający 
użytkownika w wyborze parametrów początkowych oraz kontrolujący prze- 
bieg zadania [22]. W zależności od potrzeb, kontrolowanych przez sterowanie 
KS, uruchamiany jest jeden z kilku sposobów wyznaczania kierunku poprawy 
(ryc. 9), a następnie jedna z kilku metod minimalizacji w kierunku. Punkt 
startowy zadania może być losowany lub wskazywany przez system ekspercki 
bądź bezpośrednio przez użytkownika. 


X2 


Ryc. 9. Sposoby wyznaczania kierunku poprawy w metodzie ortogonalno-diagonalnej 


Podobną formę dekompozycji można zastosować w stosunku do metody 
polioptymalizacji. Na rycinie 10 zdefiniowano sąsiedztwa S, punktu x w przy- 
padku przestrzeni N = 2 i N =3. Algorytm wyznaczający oceny niezdomino- 
wane (ryc. 11) startuje z rozwiązania narożnego ekstremalizującego funkcję 
J1(x). Bada się, które z rozwiązań należących do sąsiedztwa analizowanego 
punktu daje oceny niezdominowane. Następnie wybiera się z danego sąsiedz- 
twa rozwiązanie mające najlepszą wartość funkcji f,(x) i tworzy się wokół 
niego nowe otoczenie itd. Wystarczy założyć, że kolejne sąsiedztwo S* wektora 
XND nie musi być z góry zdeterminowane jako najliczniejsze sąsiedztwo Sy. 
Jeżeli analiza sąsiedztwa S4 dałaby w wyniku ocenę i rozwiązanie nie- 
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zdominowane, to nie byłoby konieczności przeszukiwania sąsiedztw S$, 
S5,..., SĘ. Zastosowanie takiej procedury znacznie przyspiesza proces wy- 
znaczenia rozwiązań niezdominowanych. Zdekomponowana metoda poliop- 
tymalizacji może być przy tym kontrolowana przez użytkownika lub przez 
system ekspercki, który wskaże sąsiedztwo graniczne S,, n<N, jakie ma 
powodować uaktywnienie warunku stopu. 


4;= (11,01, ZPESA URE A = [dz; Az]! 
u -1, OJ, -1, 1], 
» 3 h 1], (1, -1], 
[O, -1]) Fila 


A=[4, Az, 42] 41=([1, 0,0], A2=([1, 1,0], As=([1, 1, 1], 
[-140, ÓJ, -1, 1,0], Et. 1, AJ, 
So 5 Sr US3 (0, 1, 0), [gl+:0], [1,<1,1), 
[Ó.-<1, 0], 121,0], BGSCE 
[0-01], (1,0, 1], Mt =1t 
[0,0, -1]) (1,0, 1], [1 1421], 
IMO erij: (l,-1, —1], 
[-1, 0, -1), Fi15-1]$ 
(0, 1, 1] 
(0, -1 1] 
(0. 1, -1] 
[0 1,-1]) 


Ryc. 10. Dyskretne sąsiedztwa punktu x w zadaniu z dwoma (a) i trzema (b) zmiennymi 
decyzyjnymi 
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PREOPTE""PE. NOI ROEY. 
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f(x) 


10 20 30 40 50 t,(x) 


Ryc. 11. Wyznaczanie dyskretnych zbiorów ocen (b, d, f) i rozwiązań niezdominowanych (a, c, e) 


5 — Szczecińskie Roczniki Naukowe t. XI, z. 1 
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7. DEKOMPOZYCJA OBIEKTU OPTYMALIZACJI 


Podział konstrukcji na współpracujące ze sobą elementy, jako najbliższy 
tradycyjnemu sposobowi projektowania, jest najczęściej stosowany. W kon- 
strukcji w naturalny sposób wydzielić można elementy różniące się materiałem, 
charakterem przenoszonych obciążeń, sposobem projektowania lub techno- 
logią wykonania. Każdy z elementów konstrukcji może mieć własne kryteria 
optymalizacji f;* (x) zależne od lokalnych zmiennych decyzyjnych xf/. Wektor 
zmiennych decyzyjnych obiektu x składa się ze zmiennych decyzyjnych 
globalnych x, i lokalnych x, 


I 
CEER Gel TEdZIE+ Miz= JE sg 

i=l 
Wektory zmiennych decyzyjnych lokalnych xf? opisują indywidualne właś- 
ciwości poszczególnych elementów konstrukcji, np. rodzaje materiałów, typy 
i katalogi przekrojów, sposoby połączeń itd. Zmienne decyzyjne globalne x, 
muszą być wspólne dla co najmniej dwóch elementów, np. sposób podparcia 
konstrukcji dachowej ma bezpośredni wpływ na projektowanie przekrycia 
dachowego, słupów i fundamentów hali sportowej [3]. Obszar dopuszczalny 
zadania optymalizacji obiektu może być quasi-wielospójny lub łatwo go do 
takiego przekształcić (ryc. 7): 


I 
M =X GUZA zk" tz EX APEŻY: 
i=l 


Cechę tę wykorzystać można w algorytmie optymalizacyjnym, pozwalają- 
cym na szybsze znalezienie rozwiązania. Poszczególne podzbiory dopuszczalne 
są najczęściej różnej wymiarowości. W zależności od określenia kryteriów 
optymalizacji dla poszczególnych zadań możemy mieć do czynienia z na- 
stępującymi typami problemów zestawionymi symbolicznie w tabeli 2: 


Tabela 2. Typy problemów optymalizacji wielopoziomowej 


Zadanie Kryteria optymalizacji 
pati | |aEOB| (KA 2) W |(EZA [EÓH| 
Globalne Ja(%) Ja(x) Ja(x) fa (x) 


Lokalne JA (x) i f (x) A? (2) i (x) 
fO (x) O (x) 


Oznaczenia: f(x) — skalarna funkcja celu, f(x) — wektorowa funkcja celu 


— zadanie globalne i zadania lokalne są jednokryterialne (typ 1), 
— zadanie globalne i zadania lokalne są wielokryterialne (typ 2), 
— zadanie globalne i zadania lokalne są jedno- lub wielokryterialne 


(typ 31 4), 
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— zadanie globalne jest jednokryterialne, a zadania lokalne są jedno- 
i wielokryterialne (typ 5), 

— zadanie globalne jest wielokryterialne, a zadania lokalne są jedno- 
i wielokryterialne (typ 6). 

Zadania optymalizacji wykorzystujące dekompozycję obiektu formułowane 
są w literaturze głównie jako jednokryterialne [2, 9, 14, 17]. Jeszcze ciągle 
powstają prace dotyczące dekompozycji zagadnień liniowych, ale odniesionych 
do optymalizacji wielokryterialnej [32]. W pracy [23] podano czteropoziomo- 
wy algorytm dekompozycji w przypadku nieliniowej optymalizacji wielo- 
kryterialnej.j W artykule [28] autor przedstawia koncepcję ewoluującego 
zadania dyskretnej nieliniowej optymalizacji wielokryterialnej. 

Globalne kryteria optymalizacji formułowane są najczęściej jako minimum, 
ze względu na zmienne decyzyjne globalne, sumy minimów funkcji celu zadań 
lokalnych, określonych dla własnych zmiennych decyzyjnych lokalnych 
xf) oraz odpowiednich zmiennych decyzyjnych globalnych x, [3, 6, 17, 32]: 

I 


J(Xyp) = Min f(xg, x) = Min )) Min f(xg, xQ). 


xGEXG xGeXGi=l x(eXf) 

x(Jex() 
Jako globalne kryteria optymalizacji konstrukcji przyjmuje się najczęściej 
(tab. 1) minimum masy, objętości lub kosztu materiału, minimum kosztu 
wykonania obiektu, minimum nagromadzonej energii sprężystej, minimum 
energochłonności itd. [24]. Nie można natomiast traktować, że niezawodność 
lub sztywność całej konstrukcji jest sumą niezawodności lub sztywności jej 
elementów. Są to typowe funkcje celu zadań lokalnych, a jeżeli mają być 
przyjęte jako globalne, to muszą mieć inną postać niż sumacyjna pokazana 
wyżej [34]. 

Poszczególne zadania lokalne rozwiązywane są niezależnie, przy czym 
czynnikami koordynacyjnymi są zmienne decyzyjne globalne oraz nadrzędne 
kryteria optymalizacji. Schemat procesu optymalizacji, wykorzystujący dekom- 
pozycję obiektu opisanego funkcjami celu typu 2 (tab. 2), pokazano na ryci- 
nie 12. Podział dotyczyć może zarówno całego obiektu, jak i poszczególnych 
elementów konstrukcji. W tym drugim przypadku powstają zadania lokalne 
(L) i globalne (G) niższego poziomu. Pomiędzy zadaniami lokalnymi a global- 
nym mogą być przekazywane następujące informacje: 

— pełne zbiory ocen i rozwiązań niezdominowanych, 

— Ocena i rozwiązanie preferowane [8, 12, 28], 

— oceny i rozwiązania reprezentatywne [2, 20, 23], w przypadku licznych 
rozwiązań niezdominowanych nieznacznie różniących się ocenami. 

Pierwsze podejście, aczkolwiek najwłaściwsze z punktu widzenia popraw- 
ności założeń optymalizacji wielokryterialnej, może powodować generowanie 
zbyt wielkiej liczby globalnych zadań cząstkowych. Dla każdego lokalnego 
rozwiązania niezdominowanego musi być bowiem wyznaczony zbiór roz- 
wiązań niezdominowanych zadania globalnego. W przypadku przekazywania 
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Ryc. 12. Schemat procesu polioptymalizacji wielopoziomowej typu 2 z dekompozycją obiektu 
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wyłącznie rozwiązań preferowanych najszybciej osiągany jest wynik końcowy, 
ale kosztem utraty wielu rozwiązań niezdominowanych. Pewien kompromis 
w tym zakresie stanowi koncepcja określania i przekazywania zbioru ocen 
i rozwiązań reprezentatywnych złożonych z 2—5 elementów. 

W zadaniach polioptymalnego kształtowania złożonych obiektów mogą 
być wykorzystywane wszystkie formy dekompozycji. Poszczególne zadania 
lokalne i globalne mogą mieć zdekomponowane wektory zmiennych decy- 
zyjnych i funkcji celu. Do poszukiwania rozwiązań niezdominowanych 
można zastosować efektywne algorytmy oparte na dekompozycji metody 
optymalizacji. 

Na rycinach 13a i 13b pokazano koncepcję rozwiązania konstrukcyjnego 
lekkiej ocieplanej hali biurowo-magazynowej oferowanej przez firmę Lindab 
[18]. Zaznaczono tu wybrane globalne i lokalne geometryczne zmienne 
decyzyjne oraz wskazano możliwości podziału obiektu na niezależne elementy. 
Pełną dekompozycję hali na zadania I, II i III poziomu przedstawia schemat 
na rycinie 14. Kryteria optymalizacji, jakie można przyjąć w przypadku 
projektowania tego typu obiektów, są następujące: 

— kryterium nadrzędne: minimum oczekiwanego kosztu obiektu przy 
spełnieniu wymagań w odniesieniu do pełnionych funkcji, 

— kryteria zadaniowe: 
ekonomiczne — minimum kosztów wbudowanych materiałów, 
technologiczne — maksimum technologiczności rozwiązania systemowego, 


9.2 y2 j [z 42 2 
L= 2 *Xgy / lą *XGq  Xaz 


Ryc. 13a. Rozwiązanie konstrukcyjne hali typu Lindab 
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Ryc. 13b. Rozwiązanie konstrukcyjne hali typu Lindab 


ekologiczne — minimum energochłonności eksploatacji obiektu, 
konstrukcyjne — maksimum nośności sprężystej lub sprężysto-plastycznej, 
funkcjonalne — maksimum elastyczności funkcji obiektu. 


Poszczególne zadania lokalne mogą mieć własne kryteria optymalizacji, i to 
niekoniecznie zgodne z kryteriami zadania wyższego poziomu [28]. Przy- 
kładowo układ nośny hali, tzn. słupy i rygle ramy, można analizować według 
kryteriów: 

— minimum masy elementów konstrukcji, 
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Ryc. 14. Zadania optymalizacji I, II i III poziomu na przykładzie hali magazynowo-biurowej 


— maksimum wyrównania naprężeń w układzie o wzmocnionych węzłach, 

— minimum długości spoin (lub masy spoin, przy spoinach przerywanych 
albo o różnej grubości) niezbędnych do wykonania wzmocnienia węzłów 
narożnych, 
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izolację cieplną obiektu według kryteriów: 
— minimum kosztów materiału, sprzętu i robocizny przy wykonaniu 
ocieplenia, 
— minimum kosztów ogrzewania i wentylacji obiektu, 
a posadzkę części magazynowej według kryteriów: 
minimum kosztów materiału, wykonania i eksploatacji, 
maksimum odporności na chemikalia, Ścieranie i uderzenia, 
maksimum przeciwpoślizgowości. 


8. WNIOSKI 


1. Dekompozycja zastosowana w dyskretnym zadaniu polioptymalizacji 
konstrukcji stanowi najczęściej jedyną efektywną drogę do uzyskania zadowa- 
lającego rozwiązania. 

2. Dekompozycja, tak jak i optymalizacja, powinna być stosowana szcze- 
gólnie w rozwiązaniach systemowych i w tzw. problemach wielkich, gdzie 
nawet nieduże zyski w stosunku do rozwiązań tradycyjnych, powielane 
wielokrotnie, mogą przynieść znaczące wpływy. 

3. Bez podziału zadania liczba oszacowań funkcji celu wzrasta wykładniczo 
ze wzrostem liczby zmiennych decyzyjnych [4, 21]. Z doświadczeń autora 
wynika, że w zależności od wielkości zadania, stosując dekompozycję trzeba 
zbadać od 0,2% (dla więcej niż 10 zmiennych) do 20% (dla 2 do 3 zmiennych) 
rozwiązań należących do obszaru dopuszczalnego, aby uzyskać satysfakc- 
jonujące rezultaty. Oznacza to, że korzyści płynące z zastosowania dekompozycji 
rosną wraz ze wzrostem liczby zmiennych decyzyjnych i kryteriów optymalizacji. 

4. Nie ma uniwersalnej metody dekompozycji nadającej się do rozwiązy- 
wania dowolnego zadania z zakresu polioptymalizacji konstrukcji. Jak wynika 
z cytowanego piśmiennictwa, uwaga ta dotyczy nie tylko konstrukcji budow- 
lanych. 

5. W trakcie analizy czas obliczeń i jakość rozwiązania mogą podlegać 
kontroli zarówno przez inżyniera, jak i przez system ekspercki. 

6. Na ogół bez specjalnych procedur, skutecznych tylko w nielicznych 
przypadkach, nie ma pewności, czy wyznaczone rozwiązanie jest globalnie 
optymalne lub czy wyznaczony dyskretny zbiór rozwiązań niezdominowanych 
jest kompletny i niepoprawialny. 

7. Zastosowanie w jednym zadaniu różnych form dekompozycji znacznie 
przyspiesza uzyskanie rozwiązania. 

8. Konieczność zmiany założeń dotyczących tylko pewnych elementów 
konstrukcji jest prostsza i mniej czasochłonna w przypadku modelu zdekom- 
ponowanego niż modelu jednolitego. W tym drugim przypadku każda zmiana 
wymaga kompletnego powtórzenia analizy optymalizacyjnej całego obiektu. 

9. Poszukiwanie pełnego zbioru rozwiązań niezdominowanych jest złożone 
i trudne do zrealizowania, stąd najczęściej pomiędzy zadaniami cząstkowymi 
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przekazywane jest rozwiązanie preferowane lub 2—3 rozwiązania reprezen- 
tatywne. 

10. Dyskretnych zmiennych decyzyjnych nie powinno się zastępować 
ciągłymi, gdyż optimum zadania dyskretnego może być znacznie oddalone od 
optimum zadania ciągłego [24]. Szczególnie w przypadku większej liczby 
zmiennych decyzyjnych przybliżenie wyniku zadania ciągłego do najbliższego 
dyskretnego może okazać się bardzo czasochłonne. Już przy 5 zmiennych 
decyzyjnych do sąsiedztwa danego punktu należą 242 rozwiązania, a przy 
10 zmiennych ponad 5 tysięcy rozwiązań. W przypadku dekompozycji za- 
dania problem ten nie ujawnia się w tak jaskrawy sposób, lecz pozostaje 
aktualny. 


WITOLD MAREK PACZKOWSKI 


DECOMPOSITION OF DISCRETE POLYOPTIMIZATION TASKS 
IN BUILDING CONSTRUCTIONS 


Summary 


In the paper the problem of decomposition of a discrete polyoptimization task was presented. 
It is especially a crucial problem in the analysis of real objects with a great degree of complexity. 
Obtaining the optimal solution in such objects, with presently available technical means, depends 
to a large extent on a formulation of a task which makes use of decomposition and employing 
methods of analysis based on decomposition techniques. The following elements of an optimization 
process can be subject to decomposition: a decision variables vector, an objective function vector, 
a feasible domain an optimization method, and an optimization object. All forms of decomposition 
were discussed with special attention paid to the decomposition of an optimization object and 
function vector, as well as to optimization methods. Attention was drawn to the possibility of using 
various optimisation methods in one task, which significantly accelerates obtaining the solution. 
The concept of original algorithms of a polyoptimization method and acceptable task field was 
proposed. A hundred and twelve practicable criteria for building structures and their elements” 
optimization were presented. They were divided into six groups: superior — 9, constructional — 
35, technological — 16, economic — 22, ecological — 13, and functional — 20. The possibilities 
of task decomposition were exemplified by the office and warehouse building. 


BHTOJIBĄA MAPEK HAUKOBCKHM 


NEKOMIIO3MHIIM A 3AJAHUŃ HAACKPETHOŃ 
HOJIAONTUMMHU3ANUM CIPOHTEJIBHBIX KOHCTPYKIIUŃ 


Pe3oMe 


B pa6oTe rIper1cTaBJieHa Ipo6JeMa jKeKOMIIOZUNUH 3ANAHHA JJACKDETHOŃ HOJIMOITUMU3ALNHA. 
Ta rpoóJleMa aBJIAeTcA OCOGeHHO CyliecTBeHHOŃ HpH aHaJIu3e peAJIbHbIX OÓbEKTOB C OOJIbIIIOjń 
CTEIIEHbIo CJIOXXHOCTH. IIoJryueHHe OIITUMAJIbHOTO DpELIEeHHA TAKHX OÓbEKTOB, IIDH CylNeCTBYFO- 
LIIHX TEXHHUECKHX CPEĄCTBAX, B OOJIbIIIOŃ CTEIIEHH 3ABHCHT OT (DODMYyJIMDOBKH 3AJĄAHHA, ACIIOJIb- 
BYFOLINETO JIEKOMIIOZHIIMIO H IIDUME€HCHHE METOJIOB AHaJIH3A, OCHOBAHHbBIX Ha TEXHHKE JIEKOM- 
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IO3HIIHH. JIEKOMIIO3HUNHH MOTYT IIOJIBEpTATECA CJIERYIOINAE JJIEMEHTBI IIDOLIECCA OITHMH3ANUN: 
BEKTOD IIEpeMeHHblX, BEKTOD (DyHKIHH IIEJIM, OÓJIACTŁ JIOIIyCKAEMBIX pELIEHHi, METOĄ OII- 
THUMH3AHNAH H OÓbEKT OITAMKH3ANAA. OOCYX1EHbI Bce (©OPDM5I HIEKOMIIO3AINHH, OOÓpAaLNAA OCOÓEH- 
HOe BHHAMAHHE Ha IEKOMIIO3HIIHIO OÓbeKTA OIITHMH3AIIHM, BeKTODp QyHKNHH NEJJH H METOJ 
OITUMK3ANHH. BbIJIo OÓpalneHO BHHMMAHHe Ha BO3MOXHOCTb IIDAMCHHTb B OHHOM 3ANAaHHH 
pa3JIMUHbIie QOPMbI HEKOMIIO3HIIMM, UTO B 3HAUHTEJJbGHOi Mepe yCKODPAET IIOJIYMCHHE DELIEHHA. 
Takxe ÓblUJIM IpeHJIOXEHbI KOHLIENIIAA OPUTHHAJbHbIX AJITODHTMOB NEKOMIIO3HNHH METOJIbI 
NOJIMONTAMA3AIIHH HM OOJIACTH NOIYCKAeMbIX pelieHHi. B pa6oTe HpezyloxeHbi 112 kpuTepueB, 
BO3MOXHBIX JIJIA IIDAMCHEHHA B OITTHMH3AIIMH CTPOHTEJIbLHBIX KOHCTDYKIIMŃ MH HX 3JIEMEHTOB. OHH 
ObIJIM IIOHDPA3JIEJJEHbI HA IIECTŁ TPYylIl: TJIABHblIe — 9, KOCTPYKIIAOHHbIe — 35, TEXHOJIOTHYEC- 
kie — 16, 3KOHOMHYUECKHE — 22, 3KOJIOTHUECKHE — 13 H (KYHKIIAOHAJIbHble — 20. KblIu 
IIOKA3AHbl BO3MOXHOCTH HEKOMIIO3HNHH 3ANaHMHA Ha IIpHuMepe OQUCHO-CKJIAĄCKOTO IIOMEIIIEHHA. 
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Streszczenie: Przedstawiono jednowymiarowy nieliniowy model konsolidacji 
gruntów organicznych. Model uwzględnia zmiany parametrów gruntowych wywoła- 
ne tym procesem. W pracy omówiono również opis matematyczny badań edomet- 
rycznych gruntów organicznych oraz wskazano na sposób interpretacji uzyskanych 
wyników. Przedstawiony model może być stosowany do prognozy konsolidacji 
gruntów słabych (torfów). 


Słowa kluczowe: mechanika gruntów — konsolidacja gruntów. 


1. WPROWADZENIE 


Zagadnienie konsolidacji gruntów jest problemem, który występuje bardzo 
często w praktyce budowlanej. Podłoże gruntowe obciążone konstrukcją 
zewnętrzną podlega osiadaniu. Osiadanie gruntu odbywa się kosztem zmniej- 
szenia jego porowatości. Jeżeli pory wypełnione są wodą gruntową, to 
osiadanie — zmniejszenie porowatości gruntu — jest możliwe dopiero po 
wyciśnięciu części wody z porów gruntowych. Można zatem powiedzieć, że 
osiadanie gruntu jest opóźnione przez odpływ wody z porów. Konsolidacji 
towarzyszy pojawienie się ciśnienia wody w porach. W podłożu obciążonym, 
które podlega konsolidacji, szkielet gruntowy przenosi tylko część obciążenia 
zewnętrznego, natomiast pozostała część obciążenia zewnętrznego powoduje 
podwyższenie ciśnienia wody w porach gruntowych. 

Odpływ wody z porów gruntowych zmniejsza ciśnienie porowe. Po 
odpowiednio długim czasie można przyjąć do praktycznych obliczeń, że całe 
obciążenie zewnętrzne jest przejmowane przez szkielet gruntowy, a ciśnienie 
zmniejsza się do zera. W praktyce inżynierskiej podstawowe znaczenie ma 
znajomość czasu konsolidacji [5, 6, 14|. Zwykle przyjmuje się, iż czas ten 
odpowiada osiągnięciu 0,9 osiadania docelowego. 

Wystąpienie ciśnienia wody w porach ma istotny wpływ na warunki 
równowagi podłoża gruntowego poddanego obciążeniu. W czasie konsolidacji 
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składowe normalne naprężenia, które występują jako stan naprężenia 
w gruncie od obciążenia zewnętrznego, są pomniejszone o ciśnienie porowe. 
Fakt ten zmienia warunki równowagi granicznej w gruncie. Z tego powodu 
w praktyce inżynierskiej stosuje się wiele zabiegów, aby proces konsolidacji 
przyspieszyć poprzez szybsze odprowadzenie wody z warstwy kosolidowanej. 
W niniejszym opracowaniu analizuje się konsolidację w warunkach natural- 
nych [14]. 

Podstawowy opis procesu konsolidacji podał Terzaghi [cyt. za 14]. 
Przedstawiając opis matematyczny tego zjawiska przyjął on, że parametry 
gruntu nie zmieniają się wraz z konsolidacją. W szczególności stałe są: moduł 
ściśliwości gruntu oraz współczynnik filtracji. Badania prowadzone dla grun- 
tów organicznych [2, 3, 4, 13] wskazują, że w miarę jak postępuje konsolidacja 
i zwiększa się osiadanie gruntu, zmniejsza się jego porowatość, a wraz z nią 
zmieniają się parametry gruntu. Jednocześnie z postępującą konsolidacją 
współczynnik filtracji do wzoru Darcy maleje [2], natomiast moduł ściśliwości 
rośnie [13]. Obie te wielkości w istotny sposób wpływają na czas konsolidacji. 
Zmniejszanie się współczynnika filtracji wpływa na wydłużenie czasu kon- 
solidacji, natomiast wzrost modułu ściśliwości wpływa na skrócenie czasu 
konsolidacji. Określenie równoczesnego wpływu tych czynników na przebieg 
konsolidacji jest przedmiotem niniejszej pracy. Jeżeli badania konsolidacji 
prowadzone są w edometrze, to istotne znaczenie ma tarcie próbki gruntu 
o wewnętrzną powierzchnię cylindra edometru [8]. Tarcie to wywołane jest 
zarówno przez siły docisku, jako składowe obciążenia zewnętrznego, jak i przez 
siły adhezyjne, wynikające z fizycznych właściwości gruntów organicznych. 
Ponieważ w przypadku gruntów organicznych obciążenia pionowe nie prze- 
kraczają 100 kPa, siły tarcia mogą stanowić nawet do 20% tego obciążenia. 
W pracy przedstawiono sposób, który pozwala na uwzględnienie tego czynnika 
w interpretacji wyników badań edometrycznych. 


2. PODSTAWOWY MODEL KONSOLIDACJI GRUNTU 


Podstawowy model konsolidacji gruntu przedstawił Terzaghi [cyt. za 14]. 
Rozważania swoje przeprowadził dla jednowymiarowego układu izotropowej 
próbki gruntu zakładając, że parametry ośrodka nie zmieniają się w czasie. 
Założył dodatkowo, że pomiędzy obciążeniem a odkształceniem zachodzi 
relacja liniowa: 


! 


ds _o (1) 
dz E' 

oraz że wszystkie pory w gruncie są wypełnione wodą, a ruch wody można 
opisać jednowymiarowym równaniem filtracji Darcy (ryc. 1) 


e izernk 3x: (2) 


79 


W równaniach tych 
k, E = const. (3) 


Schematycznie analizowaną próbkę pokazano na rycinie 1. 


dz 


Ryc. 1. Schemat próbki poddanej konsolidacji 


W modelu konsolidacji Terzaghi zakładał, że wszystkie pory wypełnione są 
wodą, a zmiana porowatości wybranej objętości kontrolnej gruntu wynika 
z odpływu wody. Zgodnie z ryciną 1 mamy 


dsA = [V,(z+ dz, t)—V(z, t)] Adt. (4) 
Równanie to można również zapisać w postaci: 
O (ds ÓW, 
Aj zA |= są 5 
0i (Z) óz ©) 


Podstawiając do powyższej zależności związki (1) oraz (2) uzyskujemy pod- 
stawowe równanie różniczkowe konsolidacji. Ponieważ 


G =0qg—u, (6) 
(oj l 0| kóu 
z 6-we|- -zie%8| % 
Terzaghi w rozważaniach swoich założył, że 
E, k, 09 = const. (8) 


Otrzymał proste równanie różniczkowe konsolidacji (równanie typu dyfuzji) 
w postaci: 


0u H*ó*u 0) 
PAR ET 
gdzie wielkość T jest stałą konsolidacji 
H* 
Paz (10) 


Ek 
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Rozwiązanie równania (10) podaje się zwykle w postaci szeregów Fouriera, 
przedstawiając funkcje: 


u=uf(z, t) oraz s=s(i). 


Przy rozwiązywaniu praktycznych przypadków konsolidacji najważniejsze 
znaczenie ma czas konsolidacji, a miarą konsolidacji jest osiadanie kon- 
solidowanej warstwy gruntu w czasie. Jeżeli osiadanie końcowe (docelowe) 


S. = lim S(t), (11) 


£—>o 
to stopień konsolidacji określa się jako 


S(t 

50 _ F(g. (12) 

Na 
W praktyce inżynierskiej za czas konsolidacji uznaje się czas potrzebny 
na osiągnięcie 0,98. Oznacza to, że czas konsolidacji określa się rów- 
naniem: 


S(Te) = 0,98. (13) 


Obliczenia czasów konsolidacji wskazują, że z dostateczną dla celów praktycz- 
nych dokładnością można w modelu Terzaghiego przyjąć 


2 


Fe = WE 
Ek 


(14) 


Jak już zaznaczono, procesowi konsolidacji towarzyszy zmiana paramet- 
rów gruntowych związanych ze zmianą porowatości. Wobec tego prawidłowe 
określenie przebiegu konsolidacji wymaga uwzględnienia w obliczeniach zmia- 
ny tych parametrów. 


2.1. WPŁYW KONSOLIDACJI NA ZMIANĘ WSPÓŁCZYNNIKA FILTRACJI 
GRUNTU ORGANICZNEGO 


Jak już wcześniej wspomniano, konsolidacja gruntu związana jest ze 
zmianą porowatości. W miarę odpływania wody z porów gruntowych po- 
stępuje proces osiadania i maleje porowatość. Zmiana porowatości wpływa na 
zmniejszenie współczynnika filtracji. W literaturze [14] podaje się zależność, 
która uwzględnia zmianę współczynnika filtracji wywołaną zmianą porowato- 
Ści. Zależność ta ma postać: 


k=ce". (15) 


Szczegółową analizę zmian współczynnika filtracji wywołaną zmianą porowa- 
tości w gruntach organicznych przeprowadzili Dereczenik i Seul [2]. Po 
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podstawieniu do wzoru (15) odpowiednich zależności wyrażających zmianę 
wskaźnika porowatości, za pomocą wielkości osiadania warstwy gruntu 
poddanej konsolidacji, otrzymali związek 


S$ x 
k=[9)=kof1 5) (16) 
W równaniu (16) występują trzy parametry, które należy wyznaczyć doświad- 
czalnie dla danego rodzaju gruntu. Można to zrobić metodami statystycznymi 
(np. najmniejszych kwadratów) mając dany ciąg zmierzonych wartości £k,, S,;. 
Badania przeprowadzone dla torfów z rejonu Szczecina wskazują, że w prze- 
ciętnych warunkach a = 0,85 oraz x = 7,0. W całym zbiorze próbek wartość 
a zawierała się w przedziale (0,7--0,9), natomiast x (5,0-- 7,5). 

Stałą k, łatwo wyznaczyć z wzoru (16). Dla S$ = 0 mamyk=kśy, czyli jest to 
współczynnik filtracji próbki przed konsolidacją. Podkreślić należy, iż badania 
laboratoryjne muszą być przeprowadzone przy założeniu, że wszystkie pory 
gruntowe wypełnione są wodą. Otrzymane wyniki badań przedstawiono na 
rycinach 2—4. 


2.2. WPŁYW KONSOLIDACJI NA ZMIANĘ MODUŁU ŚCIŚLIWOŚCI 
GRUNTU ORGANICZNEGO 


Jak już wykazano, w czasie konsolidacji zmienia się moduł ściśliwości. 
Zmniejszenie porowatości gruntu w warstwie konsolidowanej wywołuje wzrost 
modułu ściśliwości. Badania zmiany modułu ściśliwości gruntu organicznego 
przeprowadzili Meyer oraz Dereczenik [13]. Przyjęto, że zmiana modułu 
ściśliwości zależna jest od zmiany wskaźnika porowatości gruntu organicz- 
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Ryc. 2. Zależność współczynnika x od porowatości 
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Ryc. 3. Zależność współczynnika k, od porowatości początkowej 
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Ryc. 4. Wykresy funkcji optymalizacyjnych 
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nego i może być oparta na wzorze (15). Po wyrażeniu wskaźnika porowatości 
za pomocą wielkości osiadania otrzymano zależność: 


S (rz 
E (0) = (1-3) A (17) 


Równość (17) można też przedstawić jako funkcję obciążenia wykorzystując 
podstawowy związek (1). Otrzymamy 


E (0) = E rz rzz | (18) 


Na podstawie powyższego równania można obliczyć, jak zmienia się osiadanie 
warstwy konsolidowanej. Ponieważ lokalny moduł ściśliwości 


do 


E(o) = Ha, (19) 


więc po wykonaniu całkowania i podstawieniu warunku początkowego, że dla 
o=(0 mamy $=0, otrzymano 


1 o 8 
s- ak 1-(1+z=gz) | (20) 


We wzorze (17) występują dwa parametry: a oraz f. Są one charakterystycz- 
ne dla danego rodzaju gruntu organicznego i należy je wyznaczyć laboratoryj- 
nie. Analiza statystyczna wyników badań laboratoryjnych torfów z rejonu 
Szczecina przeprowadzona w Katedrze Geotechniki Politechniki Szczecińskiej 
[13] wskazuje, iż istnieje możliwość określenia ekstremalnych wartości para- 
metrów a oraz 5. Ekstrema odchyłek w metodzie najmniejszych kwadratów 
pozwalają na bardzo dokładne określenie parametru 5 (ryc. 5). Bardzo dobrą 
zgodność wykazuje również aproksymacja S(o) (ryc. 6). Na rycinie 7 przed- 
stawiono zmianę modułu ściśliwości badanych próbek jako funkcji osiadania. 
Zestawienie wyników wartości ekstremalnych parametrów przedstawiono 
w tabeli 1. 

Otrzymana wcześniej zależność (20) jest nieliniowym związkiem osiada- 
nie-obciążenie w gruntach organiczonych. Równanie to można nazwać nieli- 
niowym modelem osiadania gruntów organicznych. Nieliniowość wynika 
z tego, że moduł ściśliwości w gruntach organicznych zmienia się (rośnie) wraz 
z obciążeniem. Wykorzystując wzór (20) można teraz przedstawić zmianę 
współczynnika filtracji gruntu organicznego jako funkcję obciążenia. W wyni- 
ku przekształceń otrzymano 


1 x(1—8) 
k(a) = ko(L=DEJ | (21) 
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Ryc. 5. Wykres 0 
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Ryc. 6. Aproksymowane krzywe s 
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Ryc. 7. Wykres E = f(s) oraz E = f(o) 


Tabela 1. Wyniki aproksymacji parametrów 
" E, 6 
[kPa] [kPa]? 
101,0 351-107 
168,7 1,06:107* 
191,3 1,18:107* 


3. RÓWNANIE KONSOLIDACJI GRUNTU ORGANICZNEGO 


Przedstawione na wstępie równania konsolidacji (równania od (1) do (10)) 
otrzymano przy założeniu jednoosiowego ściskania próbki o izotropowych 
właściwościach. Równania te nie mogą być stosowane dla gruntów organicz- 
nych. W gruntach organicznych zmienia się moduł ściśliwości oraz współczyn- 
nik filtracji w miarę jak postępuje konsolidacja. Aby uwzględnić te zmiany 
gruntu organicznego, do równania konsolidacji wprowadzono zależności (18) 
oraz (21). Do rozważań przyjęto taką samą zasadę jak na wstępie. Wybieramy 
w gruncie objętość kontrolną, jej kształt jest wynikiem stanu naprężenia 
w gruncie. Zakładamy, że wszystkie pory wypełnione są wodą. Zmieniamy stan 
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naprężenia, np. przez zmianę obciążenia zewnętrznego. Kształt objętości 
kontrolnej ulega zmianie w ogólnym przypadku na skutek działania składo- 
wych normalnych i stycznych naprężenia na powierzchni zewnętrznej wybranej 
objętości kontrolnej. Przejęcie tych normalnych i stycznych naprężeń przez 
szkielet gruntowy wewnątrz objętości kontrolnej wymaga zmniejszenia objęto- 
Ści porów, a co-.za tym idzie, wypłynięcia części wody z porów na zewnątrz. 
Opis zmian objętości w przypadku, gdy odkształcenia następują wzdłuż trzech 
osi, natomiast odkształcenia postaciowe w trzech płaszczyznach układu współ- 
rzędnych, wymaga wprowadzenia wielu parametrów. Ze względu na analizę 
gruntu organicznego wystąpią tu trzy różne moduły ściśliwości wzdłuż trzech 
osi oraz trzy moduły odkształcenia postaciowego. Niezależnie od tego współ- 
czynniki Poissona, które pozwalają określić odkształcenie z kierunku, w jakim 
działa siła na dwa pozostałe kierunki, są również różne w kierunku każdej 
z osi. Występowanie tak wielu różnych stałych materiałowych komplikuje 
rozważania. Dlatego w niniejszym opracowaniu ograniczono się do rozważenia 
jednoosiowej konsolidacji gruntu organicznego. Przy tych uproszczeniach 
i założeniu, że odkształcenie następuje jedynie w kierunku osi pionowej, 
równanie konsolidacji przyjmie postać: 


SJeEH|- „sg SEL 
a(żo)- Ame, SE 


Po podstawieniu zależności (18) oraz (21) otrzymamy 


0 E, k, 0 09—u)f"" "Pau 
tajkmeN Fu mi EJ PEAK ZZ IE -| s. 
Ot ( e) "m z( > E, óż-| (23) 
1+a—— 
Ło 
gdzie 
| 


Rozwiązanie tego równania polega na znalezieniu funkcji u = u(z, t) wewnątrz 
warstwy kosolidowanej, to znaczy O<z<H. 

W ogólnym przypadku funkcja ta zależy od parametrów stałych wy- 
stępujących w równaniu (23), to znaczy: a, f, x, Eg, kg, 7„. Łatwo zauważyć, 
że funkcja ta zależy również od ay, jeżeli o, = f(t, z). Oznacza to, że na proces 
konsolidacji wpływa sposób obciążania warstwy konsolidowanej oraz rozkład 
naprężeń w tej warstwie wywołany obciążeniem zewnętrznym. Rozwiązanie 
równania (23) nie jest celem niniejszego opracowania. Z równania tego można 
jednak wyciągnąć wniosek, jak zmienia się czas konsolidacji na skutek zmiany 
parametrów gruntowych. Na podstawie wzoru (10) możemy napisać, że 


sado 
AŻ TGLTGY 


(25) 
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Po podstawieniu odpowiednich zależności na zmienne E(o) i x(o) otrzymamy 
tzw. chwilowy czas konsolidacji w postaci: 


H2 1948-70-86 
Toe (i +a -E ) | (26) 
[0) 


Jeżeli oznaczymy przez” Ty chwilowy czas konsolidacji w chwili zerowej, to 
z równania (25) mamy 


(REM 

c= ; „gh 

(0) E, ko ( ) 

Po zakończeniu konsolidacji natomiast, kiedy u=0, chwilowy czas kon- 
solidacji 

M H2 05 x(B—1)-8 
I,= 1+0— 28 
|= RĘ(="E) pa: 


Z dostateczną dokładnością dla celów praktycznych można założyć, że średni 
czas konsolidacji, to znaczy czas przebiegu konsolidacji, jest średnią geomet- 
ryczną obu tych wielkości. Stąd 


gafa: 
T< mięp| PDZ Aec)|| (29) 
2 Eg 
Wzór (29) można przedstawić w postaci przybliżonej, jeżeli założyć, że 
00 
RODZA 30 
ag, < (30) 


Warunek ten występuje zazwyczaj w praktycznych obliczeniach. Otrzymamy 


wtedy 
x(B—1)—$B oc 
I ="I5,ezp 2a(8-1 | (31) 


Przekształcając ten wzór można również napisać: 


x-_l | do 
B= ly ——————- | |. 32 
: || 2a zane. G2) 
Obliczenia o zasięgu praktycznym wskazują również, że wartość parametru 
P jest zwykle duża i wynosi ok. 6,0, natomiast parametr x jest równy ok. 5,0. 


Pozwala to na dalsze uproszczenie wzoru (32). 
Ponieważ 


|| 
—1>——, 33 
K ?*pI (33) 


dla celów praktycznych można przyjąć czas konsolidacji 
K—l a 


0 
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W porównaniu z czasem konsolidacji występującym w modelu Terza- 
ghiego, otrzymany tu wynik wskazuje, iż czas ten jest zależny od obciążenia 
zewnętrznego og. Można obliczyć, że w przeciętnych warunkach, gdy o0,/E, = 
= (0,3 oraz k=50 1a=0,85, otrzymamy T = 2,42T,. 

Jak widać, uwzględnienie w obliczeniach faktu, iż w procesie konsolidacji 
parametry zmieniają się, powoduje wydłużenie czasu konsolidacji w stosunku 
do czasu obliczanego dla stałych parametrów. 


4. INTERPRETACJA BADAŃ EDOMETRYCZNYCH 
WYKONANYCH DLA TORFÓW 


Próbki gruntów umieszczane w edometrze, pobrane w terenie z pewnej 
głębokości, nie są gruntem wstępnie nieobciążonym. Zwykle próbki te poddane 
już były pewnym obciążeniom ściskającym. Jeżeli ciśnienie prekonsolidacji 
stanowi znaczącą część późniejszego obciążenia zasadniczego, to interpretacja 
wyników pomiarów edometrycznych powinna je uwzględnić. Następnym 
elementem, który powinien być uwzględniany przy interpretacji badań edomet- 
rycznych, są naprężenia, jakie pojawiają się na wewnętrznej powierzchni 
cylindra edometru. Naprężenia te przeciwdziałają obciążeniu głównemu (pio- 
nowemu), zmniejszając osiadanie próbki. Powodowane są one siłami tarcia 
wywołanymi przez poziomy docisk gruntu do wewnętrznej powierzchni 
cylindra edometru oraz przez siły adhezyjne pojawiające się na tej powierzchni. 
W próbce torfu przyjmuje się, że ma ona właściwości plastyczne oraz sprężyste. 
Założono, że odkształcenia plastyczne i sprężyste próbki o wysokości H są 
równe odpowiednio: 


H 
czaki (35) 
Ex 
s Ho 
DEM = E_ (36) 


Do dalszych obliczeń przyjęto, że założone osiadania tak w części plastycznej, 
jak i sprężystej są funkcjami liniowymi: 


HG |. 
Spl = pry => g (37) 
Oraz 
Ho ,, 
spi = E = K O, (38) 


natomiast osiadanie łączne jest Średnią ważoną: 


ć 1 


IA PERS 12] 
WIEN ANETA 


(39) 
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Z porównania zależności (35), (36) i (39) wynika, że 


y= Be, (40) 
Egpr = E(1+6). (41) 


Przyjęcie takiego założenia eliminuje w dalszej części potrzebę oddzielnego 
obliczania modułów sprężystości dla części odkształceń sprężystych i plastycz- 
nych. Metoda ta pozwala na przyjęcie jako modułu ściśliwości wielkości 
wypadkowej ze ściskania w edometrze. 

Właściwości plastyczne gruntu ujawniają się w zasadzie w sytuacji, kiedy 
odciążamy próbkę. Wówczas odkształcanie plastyczne pozostaje bez zmian, 
natomiast zmniejsza się odkształcanie sprężyste. Schematyczne osiadanie próbki 
gruntu w części sprężystej i plastycznej podczas badań edometrycznych 
pokazano na rycinie 8 [8]. Przebieg historii obciążenia jest następujący: w stanie 
naturalnym próbka obciążona jest ciśnieniem prekonsolidacji o,. W chwili t, 
próbka pobrana zostaje do laboratorium i w związku z tym odciążona do o, = O. 
W chwili t, próbka ponownie zostaje obciążona w edometrze obciążeniem o. Po 
ustaleniu się osiadania w chwili t,, odciążamy próbkę do zera, o, = O. Ponieważ 
poszczególne fazy zmiany obciążenia odbywają się wolno, w dalszej części 
założono, że każda z rozpatrywanych chwil czasowych obejmuje stan, kiedy 
osiadanie już się nie zmienia. Na podstawie ryciny 8 otrzymamy: 

w chwili ty 


S$ = Sa (t4) = Sspr(£1) = k* 04, (42) 
w chwili t, 
SAŻJZK0r. (43) 
NAM (t,) = 0, (44) 
S$= Śp: Gp (45) 
IKC R. 
| 
MEGM_s-0 MIENI swo |, 
Spl jt |” [ta KNEME 


Ryc. 8. Schemat obciążenia i odpowiadającego odkształcenia w analizowanym modelu [4] 
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w chwili t., zakładając, że 0 > Gy 


Spi(t3) = k* 03, (46) 
Sspr (£3) = K* 03, (47) 
St.) ="K"ax, (48) 
w chwili ty 
Spi (tą) = k* 03, (49) 
AWARII (50) 
S(t,) =p 0. (51) 


Na podstawie przedstawionych wyżej zależności łatwo zauważyć, że mierzona 
w laboratorium w edometrze wartość osiadania 


48, = S(t,)-S(t.) = k* pęzaik" (52) 
45; = S()-S(u) = Tizk'(a,—0,). (53) 


Elementem, który również wpływa na interpretację wyników badań labora- 
toryjnych, jest tarcie gruntu na pobocznicę cylindra edometru. W ogólnym 
przypadku siły tarcia opisane są wzorem: 


tz (z) = o(z) f+C. (54) 


Badania prowadzone w Katedrze Geotechniki Politechniki Szczecińskiej [8] 
wskazują, że w przypadku badania torfów, ze względu na małe obciążenia 
pionowe, wzór (19) można uprościć do postaci: 


1 (z) = const = c. (55) 


Schematycznie rozkład składowych pionowych naprężenia w próbce pokazano 
na rycinie 9. Na podstawie tego rysunku możemy napisać: 


o (z) = 03-17. (56) 


| 46. 6, 
1 WAPNE ||W 
INACECZRCH, 


| (z 
VLĄLĄLĘ LE dll LLlLlq<Ą 


'" o(z+dz) 
54 ! 


6(H) 


Ryc. 9. Schemat rozkładu składowych naprężenia wewnątrz próbki obciążonej w konsolidometrze 
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Wykorzystując zależności (18) oraz (19) otrzymujemy 


H 
1 ó 
AS; = g||0-ra. |tz (57) 
0 
a po wykonaniu całkowania 
: 2cH Ć 
AS; =<k (o -rzzn). (58) 
W analogiczny sposób uzyskamy 
k* 2cH 
48, = 1+ zl (03 — 0-5 6-2 | (59) 


Powyższe zależności umożliwiają wyznaczenie w czasie badań edometrycznych 
wszystkich poszukiwanych wielkości: k*, c, ć oraz o;,. 

Obliczenia parametrów przeprowadzamy metodą najmniejszych kwadra- 
tów. Wykorzystanie do wyznaczania stałych równań (58) i (59) wymaga 
zastosowania edometrów o różnych wysokościach, tak aby H, było jeszcze jedną 
zmienną niezależną. W obliczeniach wygodnie jest przyjąć najpierw, że równanie 
(58) jest równaniem warunkowym w metodzie najmniejszych kwadratów: 


Ś 0% 


ozszii A! (60) 


ASŃ = = CH, 0%] —— > [H?] A 


Określamy w ten sposób stałe E oraz c. 
Jako trzeci parametr obliczamy człon 


a BŁAŚ (61) 


Ponieważ nie znamy 6, nie możemy obliczyć wartości o,. Przyjmujemy 
wówczas 


) j H. „  6cH, 
[asf —As$]1 +. = Fi op) (62) 


jako drugie równanie warunkowe. Z zależności tej, tak jak dla równania 
prostej, obliczymy wartość ć, a następnie o,. 
Badania torfów z okolicy Szczecina dały następujące rezultaty: 


E=]135kPa, c=710kPa, 6=3,76, 0, = 22,5 kPa. 


Warto zauważyć, że uwzględnienie w obliczeniach ciśnienia prekonsolidacji 
oraz sił adhezyjnych edometru wpływa znacząco na obliczenie modułu 
Ściśliwości. Na podstawie zależności (58) mamy 


E- „ę( _20H 6 | 
= 45, i akiy pąęzz Sd|: (63) 
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Ten sam moduł obliczony dla stałych parametrów ma wartość 


H 


pihotono 
ZKE 


(64) 


5. PRZYKŁADY PRAKTYCZNEGO WYKORZYSTANIA 
PREZENTOWANEGO MODELU 


Przedstawiona w poprzednich rozdziałach metoda daje dobre wyniki przy 
obliczaniu osiadania obiektów inżynierskich (nasypy, drogi, małe obciążenia 
użytkowe) posadowionych na podłożu, w którym znajdują się grunty organicz- 
ne (torfy). Silnie nieliniowy związek pomiędzy obciążeniem a osiadaniem oraz 
duże odkształcenia warstwy słabej w gruncie sprawiają, że metoda potwierdziła 
się w praktyce. Podstawowy schemat obliczeniowy opiera się na sprawdzaniu 
warunków równowagi wybranej kolumny gruntu słabego podczas kolejnych 
faz obciążenia przy uwzględnieniu wyporu wody gruntowej. Schemat obciążeń 
kolumny przyjęty do obliczeń pokazano na rycinie 10. 


STAN 2 
STAN 1 
STAN 0 
PIERWOTNY POZIOM TERENU ĄS 
ZAYPACY, NZAÓDNY g 
| ao R 
AS 
ho 

hą ha 


POP EOEELE 


PA 


Ryc. 10. Schemat obciążeń kolumny gruntu przyjęty do obliczeń 


Podstawową relację opisującą osiadanie kolumny gruntu słabego (torfu) 
określa się zależnością (20): 


|- ady 


Zgodnie z ryciną 10 zakładamy, że kolumna gruntu w stanie nieobciążonym 
ma wysokość hy oraz moduł sprężystości E,. Następnie zakładamy, że w celu 
skonsolidowania gruntu słabego wykonujemy nasyp o miąższości H. Przy 
czym nasyp ten częściowo znajduje się w wodzie gruntowej, ponieważ na 
skutek osiadania wysokość kolumny zmniejszyła się do wielkości h,. Na 
rysunku ten stan oznaczono symbolem 1. Następnie tak skonsolidowany grunt 
słaby obciążamy dodatkowo naprężeniem Ao. Może ono oznaczać obciążenie 


93 


użytkowe lub podwyższenie nasypu. Powoduje to dalsze osiadanie kolumny do 
wysokości h,. Stan ten na rysunku ma symbol 2. Stosownie do ryciny 10 
możemy napisać następujące zależności opisujące obciążenia 1 odkształcenia 
w poszczególnych fazach: 


w stanie 1 
0, =H;),+H2(1—n)(7;—7,,): (66) 
l agiurE7 M 
Śs=h |=" = O 
ZGU=CM (r 
w stanie 2 
0, = 40+H1y, +H;(1—n)(y,—7,). (68) 


Ponieważ zakładamy, że nasyp jest nieściśliwy 


stąd otrzymamy 
0, = 01+40—48[7, —(1—n)(7,—7,)]. (70) 
Do dalszych obliczeń wygodnie jest przyjąć oznaczenie: 
), =7—(l—n)(v;—Yw)> (71) 
a zatem 
0, = 0,+40—ASy.. (72) 


Odpowiednio osiadanie wynosi 


1 Ac—ASy "3 
0,+40— 3 | (73) 


S+AS=h 1-|1— 
+ 0 a| ( + e 1 ż E, 
Odejmując stronami równania (67) oraz (73) otrzymamy związek pomiędzy 
stanem 1 oraz 2: 


1 6, 1-8 i watdcyse (i 
USZU ki | = ARŻĘARĘ | (4 
„at( +p=T R) [77 GE, | (74) 


Praktyczne wykorzystanie powyższych równań sprowadza się zwykle do 
dwóch zagadnień [11]: 

1) znamy początkową grubość warstwy słabej hy, poszukujemy natomiast 
osiadania w poszczególnych fazach konsolidacji, 

2) znamy wielkości geometryczne stanu 1 oraz 2, nie znamy natomiast 
miąższości początkowej hy. 

Dla pierwszego zagadnienia problem sprowadza się do rozwiązania rów- 
nania (74). Poszukiwana wielkość AS występuje zarówno z lewej, jak i z prawej 
strony równania (równanie uwikłane). Z tego powodu rozwiązanie równania 
wygodnie jest przeprowadzić metodą kolejnych iteracji. W pierwszej iteracji po 
prawej stronie równania (74) podstawiamy AS=0 i obliczamy osiadanie 
z pierwszej iteracji — 4AS(1). 
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Następnie po prawej stronie równania (74) podstawiamy 4ŚS = 4S(1) 
i obliczamy osiadanie z drugiej iteracji — 4AS(2). 

Ogólnie możemy napisać, że osiadanie z następnej iteracji obliczamy 
podstawiając po prawej stronie równania odpowiednią wartość z poprzedniej 
iteracji. Mamy 

AS (i) =/f[AS(i—1)]. (75) 


Okazuje się, że powyższa iteracja jest bardzo szybko zbieżna. Praktycznie 
wystarczą trzy iteracje, aby otrzymać wynik z dokładnością wystarczającą 
w obliczeniach inżynierskich. 

Poniżej podano przykład obliczeniowy osiadania torfu na Ostrowie Gra- 
bowskim w Szczecinie [11]. Do obliczeń przyjęto następujące dane: 


y,=17kN/m, y,=26kN/m, 0,=40kPa, 4o=100 kPa, 
E,=171kPa, n=023, h,=8$0m, a=0,88, 58= 146. 
Z, kolejnych iteracji otrzymano wartości 4S: 
A4S(1)=139m; 4S(2)=129m; 4A$= 1,30 m. 


Dla porównania, osiadanie bez uwzględnienia wyporu wody wynosiłoby 
1,39 m. 

W przypadku gdy rozpatrzymy zagadnienie drugie, zwykle problem for- 
mułujemy w taki sposób, aby na podstawie badań terenowych uzyskać dane, 
które pozwalają na napisanie układu równań warunkowych [9]: 


h, = f(0;, ho, a, B, E4). (76) 


Obciążenie o, dla każdego profilu badawczego opisuje się równaniem (66) lub 
(68), natomiast wysokości kolumny po zakończeniu osiadania zależnością: 


hp=h ESCE| sj 
= pa "FTaz) | (77) 


Obliczenia numeryczne wskazują, że równania warunkowe (77) pozwalają na 
ustalenia wielkości hy oraz uśrednionej dla całego obszaru wielkości Ey. 
Wykorzystujemy w obliczeniach metodę najmniejszych kwadratów odchyłek, 
przy czym określenie h, przeprowadzamy na dwóch poziomach. Do dalszych 
obliczeń wygodnie jest przyjąć podstawienie 


1-8 
-a|-(r75) |-xa Eg). (78) 


Równanie podstawowe ma zatem postać: 


h=h,X(o, Ev). (79) 


a stąd równanie warunkowe 
(80) 
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gdzie 
X, =X(o,, Eg). (81) 


Obliczenia prowadzimy w następujący sposób. Zakładamy wartość modułu 
ściśliwości gruntu słabego E, = E,. Następnie obliczamy zmienne niezależne 


X; SE X (o;, E,). 
Metodą najmniejszych kwadratów obliczamy hy: 
2 (h, X.) 
h, =h,(E,) = zXŻ * (82) 
gdzie 
X, =X(0,, E,), (83) 


oraz sumę kwadratów odchyłek 
0, =6(E) =) (h;—ho X)”. (84) 


Zmieniamy wartość modułu ściśliwości E, = E,,, 1 ponownie obliczamy h, = 
= ho(E,,,) oraz ó(E,,,). Otrzymujemy w ten sposób zbiór sumy kwadratów 
odchyłek 0, w zależności od przyjętej wartości modułu ściśliwości. Jako 
wielkość miarodajną przyjmujemy tę wartość modułu ściśliwości E,, dla której 
suma kwadratów odchyłek 0 jest minimalna. Znając miarodajną wielkość E5, 
można przejść do obliczenia pierwotnej miąższości warstwy gruntu słabego 
(torfu) zgodnie z wzorami (82) i (83). Posługując się tą metodą przeprowadzono 
obliczenia miąższości h, dla wspomnianego wcześniej Ostrowa Grabowskiego 
w Szczecinie [9]. 

W obliczeniach uwzględniono dane z 8 profilów gruntowych. Do opty- 
malizacji przyjęto następujące parametry: a = 0,87 oraz 5=15. 

W wyniku obliczeń z zastosowaniem modelu nieliniowego otrzymano 
E4=176kPa oraz h, = 11,40 m. Stosując związek liniowy pomiędzy ob- 
ciążeniem a osiadaniem (1) otrzymano E, = 213 kPa oraz h, = 11,93 m. 

Uzyskane z zastosowaniem modelu liniowego i nieliniowego wartości nie 
odbiegają istotnie od siebie. Jednak różnice wyników obliczeń mają zna- 
czący wpływ na prognozę dalszego osiadania podłoża, np. pod wpływem 
obciążenia użytkowego. Jeżeli w przypadku Ostrowa Grabowskiego założyć, 
że przykładamy dodatkowo obciążenie użytkowe Aa = 60 kPa, to licząc za 
pomocą modelu liniowego otrzymamy dodatkowe osiadanie 4S = 2,70 m, 
a za pomocą modelu nieliniowego 1,90 m. Szczegóły obliczeń przedstawiono 


w pracy [9]. 


6. PORÓWNANIE NIELINIOWEGO MODELU KONSOLIDACJI 
Z BADANIAMI INNYCH AUTORÓW 


Podstawowe znaczenie w opisie stanu gruntu ma relacja napręże- 
nie-odkształcenie. Jeżeli założyć, że pod pojęciem odkształcenia w celu uprosz- 
czenia analizy przyjmuje się względne skrócenie wysokości kolumny gruntu, 
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to związek ten zwykle przedstawia się w postaci: 


ds d 

A RE (85) 

ho, E(o) 
Model gruntu na ogół określany jest przez moduł ściśliwości E(o). Analizę 
zależności opisujących podstawowe związki E (o) przeprowadził Den Haan [3]. 
Najczęściej wykorzystywaną zależnością dla gruntów mineralnych jest wzór 


Janbu: 
o m 
E(o) = E, (z) | (86) 


gdzie E, jest założoną stałą, m — wykładnikiem potęgi oznaczonym empi- 
rycznie. 

Okazuje się, że związek ten dobrze przybliża zmiany modułu ściśliwości dla 
gruntów mineralnych, np. piasków pod dużymi obciążeniami. 

Den Haan [3] na podstawie własnych badań gruntów organicznych 
proponuje zależność opisującą związek obciążenie-osiadanie w postaci: 


8 o — const a (87) 
ah, | const | 


przy czym zastrzega się, że stała, która ma miano ciśnienia w tym wzorze, ma 
zawsze wartość dodatnią. Zależność ta opisuje dobrze związek między wskażź- 
nikiem porowatości a obciążeniem jedynie dla dużych wartości o. Dla małych 
wartości o licznik ułamka po prawej stronie równania jest mniejszy od zera 
1 w konsekwencji rozwiązanie prowadzi do liczb urojonych i jest nierealne. 

Podstawowa zależność S = S(o) dla torfów była przedmiotem szerokiej dys- 
kusji na międzynarodowym sympozjum naukowym zorganizowanym w Delft 
w Holandii w 1993 r. (Workshop: „Advances in Understanding and Modelling 
the Mechanical Behaviour of Peat”). Na sympozjum tym przedstawiono trzy 
referaty z Katedry Geotechniki Politechniki Szczecińskiej. Prezentowano w nich 
nieliniowy model konsolidacji omawiany w niniejszym opracowaniu. Analiza 
modelu przedstawionego przez Meyera [3] może być sprowadzona do 
konwencji proponowanej przez Den Haan; otrzymamy wówczas zależność 


1 o + const A (88) 
ahg ( const | 


m=|1—$5 oraz const = (8—1)aE,. (89) 


gdzie 


W powyższym równaniu stała jest zawsze wielkością dodatnią i podobnie jak 
we wzorze Den Haan ma wymiar ciśnienia. Dodatnia wartość tej stałej po- 
woduje, że zależność obciążenie-osiadanie jest ważna w całym zakresie obciążeń. 
Weryfikacja tej zależności przeprowadzona na podstawie badań eksperymental- 
nych wykazuje bardzo wysoką zgodność. Analiza sposobu wyprowadzania 
wzoru przez Den Haan oraz Meyera wskazuje, iż uzyskanie sumy składników 
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stałej we wzorze (89) wynika z uwzględnienia dodatkowych uwarunkowań 
próbki. Próba umieszczona w edometrze zwykle podlegała już jakimś ob- 
ciążeniom i w stosunku do gruntu absolutnie nie obciążanego w edometrze 
(virgin peat) mierzymy jedynie różnicę — nadwyżkę osiadania AS ponad już 
wcześniej istniejące odkształcenie S,. 

Wskaźnik porowatości próbki w stosunku do zmierzonego w edometrze 
osiadania wyraża się jako 


S 
No = 
H 
Ś) = , 90 
ŚĆ) TEG (90) 
Ponieważ w edometrze mierzymy jedynie różnicę osiadania 
otrzymamy 
( z) AS 
s |= 
H) H 
e(5) = (92) 
1-n, 
Jeżeli przyjąć podstawienie przyjęte w niniejszej pracy: 
So 
=Nno——, 93 
a = No— 77 (93) 
to ostatecznie otrzymamy 
a AS 
S)=-——11-—|. 94 
e(5) xl A (94) 


Kolejne relacje stanowiące zasadnicze związki w niniejszej pracy wynikają 
z przekształceń oraz scałkowania podstawowej zależności obciążenie-osiada- 
nie. W swojej pracy doktorskiej z 1995 r. Den Haan przyznaje, że proponowa- 
na zależność rozwiązuje problem w całym zakresie obciążeń oraz że cytowane 
u niego wcześniej ciśnienie może przyjmować wartość ujemną. Zastrzega, że 
wyciągnięcie ostatecznych wniosków wymaga dalszych badań. 


7. WNIOSKI 


Przedstawiono nieliniowy model konsolidacji gruntów organicznych. Mo- 
del ten został opracowany w Katedrze Geotechniki Politechniki Szczecińskiej 
pod kierunkiem Z. Meyera. W porównaniu z modelami prezentowanymi do tej 
pory w piśmiennictwie, omawiany model może być zastosowany w całym 
zakresie obciążeń. Zweryfikowano go w badaniach laboratoryjnych i tereno- 
wych, a uzyskane wyniki wskazują na bardzo dobrą zgodność osiadania 
obliczonego i zmierzonego. 

Program dalszych badań przewiduje sprawdzenie przedstawionego modelu 
dla gruntów torfowych o innym wieku 1 innej strukturze w celu ustalenia, czy 


7 — Szczecińskie Roczniki Naukowe t XL, z. I 
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proponowana ogólna zależność obowiązuje w całym zakresie obciążeń oraz czy 
można stałe tego modelu uzależnić od parametrów fizycznych gruntu. Celem 
badań jest opracowanie programu obliczeniowego, który dla większego obszaru 
torfów obciążonych pozwoli na prognozę osiadań w zależności od obciążeń. 


8. SPIS OZNACZEŃ 


a — parametr równania zmian modułu ściśliwości współczynnika filtracji 
e — wskaźnik porowatości 

E — moduł ściśliwości 

E, — moduł ściśliwości gruntu bez obciążenia 

H — wysokość próbki, grubość warstwy konsolidowanej 
k — współczynnik filtracji gruntu 

k, — współczynnik filtracji gruntu bez obciążeń 

k* — stała podatności podłoża 

n — porowatość 

m — wykładnik potęgi we wzorze Janbu 

$  — osiadanie 

t  — czas 

FT — stała czasowa konsolidacji 

TI, — czas konsolidacji 

IT, — chwilowy czas konsolidacji w chwili początkowej 
I, — chwilowy czas konsolidacji 

u — ciśnienie wody w porach 

+  — prędkość przepływu wody w kierunku osi z 

z — pionowa oś współrzędnych 

BP — parametr równania zmian modułu ściśliwości 

y+ — cięzar objętościowy gruntu masy 

y;, — ciężar właściwy gruntu nasypowego 

yw — ciężar właściwy wody 

k — parametr równania zmian współczynnika filtracji 
o, — składowa normalna naprężenia 

o — naprężenia efektywne 

0, — obciążenia przyłożone zewnętrznie 

1 — składowa styczna naprężenia 


ZYGMUNT MEYER 


A MODEL OF ORGANIC SOIL CONSOLIDATION 


Summary 


The development of industry and lack of grounds for new buildings forces us more and more 
often to use areas of poor soils. Poor soils, mainly peats and alluvial deposits, occur mostly in estuaries 
in river valleys. They are attractive sites for industrial buildings, because they enable the unrestrained 


99 


access to a river as a waterway and, at the same time, taking advantage of water transport. Poor 
soils are characterised by high heterogeneity of the medium and high setllement. One of the 
methods of improving these soils” characteristics is their consolidation. Consolidation is done by 
overloading a layer of a poor soil with a filling of silt. 

The model of organic soil consolidation during changes in the soil parameters in the course of 
consolidation was presented in the paper. The carried out research shows that in the process of 
consolidation the modulus of compressibility and the filtration coefficient change as the result of 
decrease in porosity of the medium. These are the two soil parameters that decide about the time of 
the consolidation. A mathematical description of these parameters and generalised constitutive 
equation of organic soil consolidation under these conditions were presented in the paper. It was 
stated on the basis of the solutions” analysis of that equation how much is the time of consolidation 
in those conditions longer in comparison to results of calculations made by classical methods. 
Apart from that, it was presented how a target soil medium settlement will change under an 
external load if the "medium strengthening” presented in the proposed model is taken into 
consideration. 


3ZBITMYHT MEEP 
MOJEJIE KOHCOJIMAANHMH OTPAHMUEHHDBIX IPYHTOB 


Pe3roMe 


Pa3BHTHE IPOMBIIIIJIEHHOCTH, A TAKXXE OTCYTCTBHE TEppUTOpNHA JUIA HOBbIX CTPOHTEJIBHBIX 
OGŁEKTOB 3ACTABJIAIOT HAC BCE HALIIE KCIIOJIL3ZOBATb IIJIOLNAJA CO CJIAÓBIMU TDYHTAMU. CJiaObie 
TDYHTBI, TJIABHbBIM OOÓpa30M TOp(bbI M HAaHOCHBIIŃ MJI, PACIIOJIOXKCHbI UALIE BCETO Y YCTbA peK 
B DpeUHbIX JIOJIMHaX. DTO IIpUBJIEKATEJIBHbIE MECTA „JIA IIDOMBIIIIJIEHHOTO CTPOHTEJIbCTBA, T.K. 
NAIOT BO3MOXHOCTb CBOOÓONHOTO JNOCTylia K peKe, B KAHUECTBE BOJNHOTO IIYTH AH TEM CAMbIM 
ACIIOJIBBOBATb BOJIHbIi TpaHCIIOpT. CJIaObie TpyHTbi XapakTepH3YFoTCA OOJIbILIOŃ HeONHOPOA1- 
HOCTbIo pańoHa HM OOJIBIIIO4 BOCIIDHMMUKBOCTBIO K HaTpy3kaM. ONUH U3 METONOB YyJIYUNIEHHA 
CBOŃCTB 3TUX TDYHTOB — 3TO KOHCOJIMNAIIAA. KOHCOJUNAIIAA ZAKIIIOUAETCA B IIeperpy3ke CJIOA 
cjiaOOoro TpyHTa IIyTeM HaCBlIIi. 

B pa6oTe npe4cTaBJIeHa MOJIEJIb KOHCOJIA NAIIMM OTDAHHUCHHBIX TDYHTOB B YCJIOBHAX, KOTA 
napaMeTpbi TDyHTa H3MEH4AIOTCA IIO Mepe KaK HaCTyllaeT KOHCOJIANAIIMA. IIDpOBENEHHbie HC- 
CJIETOBAHHA IIOKA3BIBAIOT, UTO B IIpOLIECCE KOHCOJIM HAIIMK, B CJIEJJCTBAE YMEHBIIIEHHA IIODHCTOCTA. 
pafoHa, U3MeH4eTcA MOJNYJIb CKHUMAEMOCTH TpyHTa u KkoaphunueHT (UJIBTparuM — N4Ba 
rapaMeTpa TpyHTa, KOTODpble peliiaroT O BpeMeHH KOHCOJIMIAIIAM. B pa6oTe npe1cTaBJIeHO 
MaTeMaTHUeck0e OIIMCAHHe HU3MeHCHHŃ 3THUX IIApaMETDOB M OGOÓLIEHO KOHCTUTYEHTHOE ypaB- 
HeHHe KOHCOJIM/IAIHMH OTDAHHUCHHbIX TDYHTOB B 3THX YyCJIOBHAX. Ha OCHOBe AaHaJIM3a pELIEHHA 
3TOTO ypaBHEHHA YyCTAHOBJIEHO Ha CKOJIEKO YBEJIMUMUBAETCA BPEMA KOHCOJHĄANAM B 3THX 
YCJIOBHAX IIO CDABHCHUKO C pEZYJIbGTATAMM BbIUKCJIEHMŃ KJIACCHUECKAMH METOJNAMH. He3aBACHMO 
OT 3TOTO ObIJIO OIIDENEJJEHO KAK H3MCHACTCA B KOHEUHOM MTOTE OCEJTAHHE TPYHTA IO BHEIIHCŃ 
Harpy3kKoń, €CJIM YHMTbIBATb „„,yKpEIJIEHHE pańoHa ”, IPÓRCTABJIEHHOTO B IIpeJJIOKEHHOŃ MOJEJIM. 
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BOGDAN AMBROŻEK 


MODELOWANIE PROCESÓW ADSORPCJI I DESORPCJI 
MIESZANINY DWUSKŁADNIKOWEJ 
Z NIERUCHOMEGO ZŁOŻA ADSORBENTU 


Instytut Inżynierii Chemicznej i Chemii Fizycznej Politechniki Szczecińskiej, 
al. Piastów 42, 71-065 Szczecin 


Streszczenie: Przedstawiono wyniki analizy teoretycznej procesów adsorpcji 
i desorpcji mieszaniny par n-butanolu i tetrachlorometanu w nieruchomym złożu 
węgla aktywnego. W obliczeniach zakładano, że adsorpcja prowadzona jest ze 
strumienia powietrza, natomiast regeneracja złoża węgla aktywnego strumieniem 
ogrzanego azotu w obiegu otwartym. Badania wykonano za pomocą modelu 
równowagowego, uwzględniającego wymianę ciepła przez ściankę kolumny adsorp- 
cyjnej. Równowagę adsorpcji opisano rozszerzonym równaniem Langmuira. Badano 
wpływ stężenia par rozpuszczalników w powietrzu na wlocie do kolumny na 
przebieg procesu adsorpcji oraz wpływ temperatury azotu na wlocie do złoża, 
średnicy kolumny adsorpcyjnej oraz intensywności strat ciepła przez Ściankę 
kolumny na przebieg desorpcji. 


Słowa kluczowe: adsorpcja w układzie dwuskładnikowym — desorpcja 
w układzie dwuskładnikowym — modelowanie matematyczne adsorpcji — modelo- 
wanie matematyczne desorpcji. 


1. WPROWADZENIE 


W wielu gałęziach przemysłu zużywane są duże ilości rozpuszczalników 
organicznych, które w wyniku różnych procesów produkcyjnych odparowywa- 
ne są do otaczającego powietrza. Duże znaczenie ma więc oczyszczanie 
powietrza odlotowego i odzyskiwanie z niego rozpuszczalników. Do oczysz- 
czania powietrza zanieczyszczonego parami rozpuszczalników stosowane są 
najczęściej następujące metody: absorpcja, adsorpcja, kondensacja, spalanie 
termiczne i katalityczne, metody membranowe [22]. 

Wśród wymienionych duże znaczenie mają metody adsorpcyjne. Do 
usuwania rozpuszczalników stosuje się jako adsorbent węgiel aktywny, chociaż 
znane są ostatnio przykłady zastosowania do tego celu również sit molekular- 
nych [23]. Adsorpcja prowadzona jest w adsorberach ze złożem stacjonarnym, 
obrotowym lub fluidalnym [ 14, 22]. Najczęściej używane są jednak adsorbery 
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ze złożem stacjonarnym. Instalacja adsorpcyjna zawiera zwykle dwa lub 
większą liczbę adsorberów. Pełny cykl pracy instalacji składa się z etapów 
adsorpcji, desorpcji, suszenia (w przypadku, gdy do desorpcji stosowana jest 
para wodna) i chłodzenia złoża. Zanieczyszczone powietrze przepuszczane jest 
przez złoże węgla aktywnego w jednej z kolumn. Adsorpcja prowadzona jest 
najczęściej do momentu uzyskania na wylocie z kolumny założonej wartości 
stężenia usuwanych rozpuszczalników. Równolegle z adsorpcją realizowana 
jest regeneracja węgla w drugim adsorberze. Najpierw przeprowadza się 
desorpcję zaadsorbowanych na złożu węgla rozpuszczalników. Desorpcja 
odbywa się za pomocą pary wodnej lub ogrzanego gazu obojętnego. Po 
zakończeniu desorpcji złoże chłodzi się. Rozpuszczalnik odzyskiwany jest 
w postaci skroplin po ochłodzeniu produktów desorpcji. 

Przy zastosowaniu do desorpcji pary wodnej złoże węgla w kolumnie 
powinno być osuszone ogrzanym gazem, np. powietrzem. Jeżeli desorbowany 
rozpuszczalnik rozpuszcza się w wodzie, produkty desorpcji muszą być 
poddane rozdziałowi metodą destylacji lub rektyfikacji. Poważnym problemem 
jest również konieczność odprowadzania Ścieków zanieczyszczonych rozpusz- 
czalnikami. 

W wyniku desorpcji gazem obojętnym uzyskuje się czysty rozpuszczalnik 
lub mieszaninę rozpuszczalników. Desorpcja może być w tym przypadku 
prowadzona w obiegu otwartym lub zamkniętym. Przy desorpcji w obiegu 
zamkniętym gaz obojętny przepływa przez podgrzewacz i kolumnę adsorpcyj- 
ną. Gaz po desorpcji kierowany jest do skraplacza, gdzie zdesorbowany 
rozpuszczalnik częściowo wykrapla się i następnie za pomocą wentylatora jest 
kierowany ponownie do podgrzewacza i kolumny adsorpcyjnej [4, 20]. Znane 
jest również rozwiązanie, w którym w czasie desorpcji rozpuszczalnik oddziela 
się od gazu obiegowego za pomocą membran (proces CAMP [22]). W czasie 
desorpcji w obiegu otwartym gaz desorbujący przepływa przez podgrzewacz 
i kolumnę adsorpcyjną. Gaz po desorpcji przechodzi przez skraplacz, gdzie 
rozpuszczalnik częściowo wykrapla się i następnie miesza z oczyszczanym 
gazem kierowanym do drugiego adsorbera [16]. 

Znane są również metody polegające na połączeniu adsorpcji ze spalaniem 
katalitycznym lub płomieniowym [11]. Najpierw powietrze oczyszczane jest 
z par rozpuszczalnika w instalacji adsorpcyjnej. Po zakończeniu adsorpcji 
przeprowadza się desorpcję zaadsorbowanego wcześniej rozpuszczalnika 
ogrzanym gazem obojętnym. Zatężone w wyniku desorpcji pary rozpuszczal- 
nika poddawane są spalaniu. Dzięki takiemu sposobowi prowadzenia procesu, 
zamiast dużych objętości powietrza o małym stężeniu rozpuszczalnika, do 
spalania kieruje się znacznie mniejsze objętości o wysokim stężeniu. Ponadto 
część rozpuszczalnika można odzyskać w postaci skroplin po ochłodzeniu gazu 
po desorpcji. 

Ustalenie warunków pracy instalacji adsorpcyjnej wymaga określenia 
zależności pomiędzy stężeniem adsorbatów na wylocie z kolumny w czasie 
poszczególnych etapów pełnego cyklu i parametrami procesu. Informacje ta- 
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kie mogą być uzyskane za pomocą modelu matematycznego kolumny adsorp- 
cyjnej. 

Opublikowane dotychczas w piśmiennictwie prace poświęcone modelowa- 
niu procesów adsorpcyjnego odzyskiwania rozpuszczalników skupiają się 
najczęściej na analizie etapów adsorpcji lub desorpcji pojedynczego składnika, 
podczas gdy przemysłowe gazy odlotowe zawierają często mieszaniny rozpusz- 
czalników. 

W pracy niniejszej przedstawiono wyniki analizy teoretycznej procesów 
adsorpcji i desorpcji mieszaniny dwuskładnikowej n-butanol-tetrachlorometan. 
Zakładano, że adsorpcja prowadzona jest ze strumienia powietrza na nierucho- 
mym złożu węgla aktywnego. Desorpcję prowadzono za pomocą gazu obojęt- 
nego w obiegu otwartym. Jako adsorbent wybrano krajowy węgiel aktywny 
typu N. 


2. MODEL MATEMATYCZNY PROCESU 


Korzystając z modelu matematycznego adsorbera można wyznaczyć prze- 
strzenno-czasowy rozkład adsorbatu i temperatury w złożu. Pozwala to na 
obliczenie podstawowych wielkości charakteryzujących proces odzyskiwania 
rozpuszczalników: czasu trwania adsorpcji/desorpcji, zużycia gazu do desorp- 
cji, stężenia adsorbatów w fazie gazowej na wylocie z kolumny itp. 

Opublikowane dotychczas prace poświęcone modelowaniu procesów ad- 
sorpcji/desorpcji są nieliczne, w szczególności dotyczy to desorpcji w układach 
wieloskładnikowych. Większość z nich oparta jest na teorii równowagowej, 
w której pomijane są opory przenoszenia ciepła i masy. Najważniejsze z nich to 
prace Rhee i Amundsona [18], Pana i Basmadjiana [17] oraz Jacoba 
i Tondeura [13]. 

Dokładniejsza analiza adsorpcji wymaga uwzględnienia oporów przenosze- 
nia masy i ciepła. Znaczące wyniki badań opartych na tym podejściu 
opublikowane są w pracach Chi i Wasana [6], Garga i Ruthvena [8], 
Basmadjiana [5], Tsai, Wanga i Yanga [21], Kumara i Dissingera [15] oraz 
Schorka i Faira [19]. 

Modelowaniem desorpcji nieizotermicznej dla układu dwuskładnikowego 
etan-propan-węgiel aktywny zajmowali się Huang i Fair [9, 10]. 

Davis, McAvoy i LeVan [7] modelowali cykliczną pracę adiabatycznej 
kolumny adsorpcyjnej. 

Proces adsorpcyjnego rozdziału mieszaniny gazowej dla przypadku, gdy 
adsorbuje się jeden składnik modelowali Hwang, Jun i Lee [12]. Model 
uwzględniał opory przepływu ciepła i masy oraz zmianę prędkości przepływu 
gazu nośnego. 

Stosowany w niniejszej pracy model zakłada istnienie lokalnej równowagi 
adsorpcyjnej pomiędzy fazą gazową a adsorbentem w dowolnym punkcie 
złoża. Takie przybliżenie opisywane jest powszechnie w literaturze. W przypad- 
ku adsorberów przemysłowych o dużych wysokościach złoża opory ruchu 
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ciepła i masy wpływają jedynie w niewielkim stopniu na kształt krzywych 
wyjścia stężenia i temperatury [7]. Dodatkowo przyjęto następujące założenia: 

— można pominąć akumulację ciepła 1 masy w fazie gazowej, która jest 
nieduża w porównaniu z akumulacją w fazie stacjonarnej, 

— nie występuje dyspersja osiowa 1 promieniowa ciepła i masy, 

— nie ma 'spadku ciśnienia w złożu. 

Przy takich założeniach równania bilansu ciepła i masy dla równowagowej 
adsorpcji/desorpcji dwóch składników mają następującą postać: 


OY, 0a, | 

+0, zi = =1,2 

Gz +0, 3 (0; F="12, (1) 
POBE ROT ać 0a 4k (T— T,) 

GC,3 FO Cs zę + 3 (e, AH.) = AEG: — (2) 


Zależność pomiędzy wielkością adsorpcji i stężeniem adsorbatu w fazie 
gazowej oraz temperaturą oblicza się z równania równowagi, które można 
ogólnie przedstawić następująco: 


d4=Qa;(7,, ,T), d=G(1, V, T). (3) 


Zastosowany w niniejszej pracy model matematyczny został wykorzystany 
wcześniej do przeprowadzenia analizy desorpcji w układzie jednoskładniko- 
wym n-butanol-węgiel aktywny N oraz tetrachlorometan-węgiel aktywny N. 
Wykonano weryfikację wyników obliczeń na podstawie wyznaczonych do- 
świadczalnie krzywych wyjścia stężenia. Dla obu badanych rozpuszczalników 
krzywe uzyskane za pomocą modelu poprawnie oddawały jakościowy i iloś- 
ciowy przebieg desorpcji [2, 3]. 

Pełny cykl pracy instalacji adsorpcyjnej składa się z trzech kolejno po sobie 
następujących etapów: adsorpcji, desorpcji i chłodzenia złoża. Rozkłady 
stężenia adsorbatów i temperatury po zakończeniu każdego etapu określają 
warunki początkowe dla etapu bezpośrednio po nim następującego. W niniej- 
szej pracy analizowano przebieg procesu adsorpcji i desorpcji w pierwszym 
cyklu adsorpcyjnym. Zakładano, że proces adsorpcji prowadzony jest na złożu 
węgla aktywnego nie zawierającego badanych rozpuszczalników. 

Warunki graniczne dla procesu adsorpcji mają zatem postać: 

— warunki początkowe 


dla ©="0RLX 20; Ek, =>) 50, 
Y =T(x)=0, (4) 
T = T(x) = const, 

— warunki brzegowe 

dla. £>01%=0; T = T,, = Const, 
Y, = Val = COnst, . (5) 
Y, = V.a2 = CONSst. 
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Warunki graniczne dla procesu desorpcji można zapisać w następujący 


sposób: 
— warunki początkowe 
dlat=0i1x>0: kj = X,.(36); 
2 = Y (x), (6) 
=") 
— warunki brzegowe 
dlat>0ix=0: F= 1, = const, 
4 = di = 0, (7) 
Y, = hd = 0 
Założono stałą temperaturę fazy gazowej na wlocie do złoża. 
Po wprowadzeniu zmiennych bezwymiarowych 
x 
=: 8 
z=r (8) 
Gt 
RE YĆ (9) 
00L 
układ równań (1)—(2) można zapisać: 
GY, 0a, 
i — 0, (10) 
Ó0z 00 ÓT 
"OE 046,09 „b 0a, 0a, 
C,—+7—-—4—|(4H,—+4H,— | = —k(1T—T,), 11 
Eczę zp; OT | ld 2 8 ( ) Sr 
gdzie 
4kL 
k=—. 12 
DG (12) 


Równania (10) i (11) stanowią układ równań różniczkowych cząstkowych. 
Do jego rozwiązania zastosowano metodę linii. Złoże adsorbentu dzielone jest 
na n warstw o jednakowej wysokości. Wysokość jednej warstwy wynosi 


L,=—, (13) 


a odpowiadająca jej wysokość bezwymiarowa 


4 
a 


SĘ. 14 
= (14) 


Dla warstwy i równania (10) i (11) można zapisać w następującej postaci: 


da 8gh- ij — Tj 


ij _ 


, 5 
dz 0, Az (13) 
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dT, = T;- wl; n d i d i 
w dO (ran ge)-a-co]| 06 
[o 


b=slbysx= Ri" TERA 


Zastosowana metoda sprowadza problem rozwiązania układu trzech rów- 
nań różniczkowych cząstkowych do rozwiązania układu 3n równań różnicz- 
kowych zwyczajnych. Do rozwiązania układu równań (15)—(16) zastosowano 
metodę Rungego-Kutty czwartego rzędu z automatycznym doborem kroku 
całkowania. W obliczeniach przyjmowano n = 50. 


3. WYNIKI OBLICZEŃ 


Wykorzystując przedstawiony wyżej model równowagowy wykonano ana- 
lizę teoretyczną dwóch najważniejszych etapów pracy instalacji adsorpcyjnej, 
tzn. adsorpcji mieszaniny par n-butanolu i tetrachlorometanu ze strumienia 
powietrza na nieruchomym złożu węgla aktywnego oraz desorpcji tych 
rozpuszczalników ogrzanym azotem w obiegu otwartym. Najważniejsze włas- 
ności badanego układu zestawiono w tabeli 1. 


Tabela 1. Właściwości badanego układu 


Gęstość nasypowa adsorbentu, g, 420 kg/m? 
Ciepło właściwe adsorbentu, C, 1,29 kJ/kg:K 


Ciepło adsorpcji n-butanolu, AH; —8,89 - 10% kJ/kmol 
Ciepło adsorpcji tetrachlorometanu, AH, —5,99 - 10% kJ/kmol 
Wysokość złoża, L 0,8 m 
Temperatura otoczenia, 1, 293 K 


| Gęstość strumienia molowego gazu inertnego, G 8,892 - 107% kmol/(m*: s) 


Równowagę adsorpcji dla układu n-butanol-tetrachlorometan na węglu 
aktywnym N opisano za pomocą rozszerzonego równania Langmuira: 
Am; K;Dj 


p | 0 AMA 2 17 
e 1 +k; p, + ka P> (7 


Założono, że indeks 1 odnosi się do n-butanolu, natomiast 2 do tetra- 
chlorometanu. Zależność współczynników powyższego równania od tempera- 
tury opisują następujące wyrażenia: 

a; k;, = exp(k4+k,„/T+ GAR), 


(18) 
k; = exp (ky + k;s/T+k;q/T7). 


Wartości współczynników liczbowych w powyższych zależnościach wyno- 
szą [1]: 
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kią = 17,569, k,, = —2,019- 10%, k; = 4,195-105, k,, = 27,322, 
kis = —2,505-10*, k,, = 4,881-105, k,, = —14,094, 
ką = 1,709: 103, k,; = 3,237:105, k,, = —6,665, 

ks = —2,428-103, k,, = 9,449 - 105. 


Analizę adsorpcji wykonano przy założeniu, że początkowa temperatura 
złoża równa jest temperaturze oczyszczanego powietrza i wynosi 293 K. 
Obliczenia wykonano dla różnych kombinacji stężeń par rozpuszczalników 
w powietrzu na wlocie do adsorbera. 

Przykładowe wyniki wykonanych obliczeń przedstawiono na rycinach 
1 —3. Krzywe wyjścia stężenia i temperatury dla różnych stężeń rozpuszczal- 
ników w oczyszczanym powietrzu pokazano na rycinach 1 i 2. Na krzywej 
wyjścia temperatury widoczny jest wyraźny wzrost temperatury powietrza na 
wylocie z adsorbera. W odróżnieniu od adsorpcji pojedynczego rozpuszczal- 
nika widoczne są dwa obszary plateau. Analiza krzywych wyjścia stężenia 
wskazuje, że składnikiem szybciej pojawiającym się na wylocie ze złoża jest 
tetrachlorometan. Rozpuszczalnik ten jest wypierany ze złoża przez lepiej 
adsorbujący się składnik, tzn. n-butanol. W wyniku tego na krzywej wyjścia 
tetrachlorometanu pojawia się obszar plateau, w którym stężenie rozpuszczal- 
nika przewyższa wartość na wlocie i jest zależne od stężeń obu składników 
w oczyszczonym powietrzu. 

Rozkłady stężenia adsorbatów w momencie przebicia złoża tetrachlorome- 
tanem przedstawiono na rycinie 3. Za stężenie przebicia przyjęto stężenie 
tetrachlorometanu na wylocie z adsorbera równe 5% jego stężenia na wlocie. 

Analizę desorpcji wykonano zakładając, że kierunek przepływu gazu jest 
przeciwny niż w czasie adsorpcji. Za początkowy rozkład stężenia rozpuszczal- 
ników i temperatury w złożu przyjmowano rozkład tych wielkości w procesie 
adsorpcji w momencie przebicia złoża tetrachlorometanem. Zakładano, że 
desorpcja prowadzona jest przy użyciu czystego azotu. 

Przykładowe krzywe wyjścia stężenia i temperatury w procesie desorpcji 
przedstawiono na rycinie 4. Na krzywych wyjścia n-butanolu widoczne jest 
podwójne maksimum stężenia. W odróżnieniu od krzywych wyjścia dla 
desorpcji pojedynczego składnika [2] nie obserwuje się wyraźnie rozwiniętego 
obszaru plateau. 

Rozkłady stężenia adsorbatów i temperatury w złożu dla różnych momen- 
tów czasowych dla temperatury azotu na wlocie do złoża Ty =413K 
przedstawiono na rycinie 5. Na rycinie 6 przedstawiono krzywe wyjścia 
temperatury i stężenia par adsorbatów uzyskane dla różnych wartości tem- 
peratury azotu na wlocie do złoża przy identycznych początkowych roz- 
kładach stężenia i temperatury w złożu. Obliczenia wykonano dla następują- 
cych wartości temperatury na wlocie do złoża: T, = 343, 393 i 413K. 
Widoczny jest znaczny wpływ temperatury zarówno na poziom maksymal- 
nych wartości stężenia rozpuszczalników, jak i na kształt krzywych wyjścia. 
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Ryc. 1. Krzywe wyjścia stężenia n-butanolu i tetrachlorometanu w procesie adsorpcji, 
k=1J/(m*:'s:K, D=1m 
A — Yi = 1,62:107* kmol/kmol,, Y,, = 7,82: 107 * kmol/kmol,, 
B — Y,1 = 3.25:10* kmol/kmol,, Y,,> = 1,56:10 * kmol/kmol,, 
C — Y = 3,23:107* kmol/kmol,, Y,.2 = 1,56:107* kmol/kmol, 
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Ryc. 2. Krzywa wyjścia temperatury w procesie adsorpcji, k = 1 J/(m*:'s:K, D=1m 
Ya = 1,62:107* kmol/kmol,, X, = 7,82: 10 * kmol/kmol, 
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Ryc. 3. Końcowe rozkłady stężenia adsorbatów w złożu w procesie adsorpcji, k = 1 J/(m?:s:K), 
D=l1m 
Ya = 1,62:107* kmol/kmol,, Y,,> = 7,82:107* kmol/kmol, 
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Ryc. 4. Krzywe wyjścia stężenia (A) i temperatury (B) w procesie desorpcji, k = 1 J/(m*:s:K), 
D=lm, 1,=413K 
Ya1 = 3,25:10 kmol/kmol,, Y,,2 = 1,56:10* kmol/kmol, 
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Ryc. 5. Krzywe rozkładu stężenia n-butanolu (A), tetrachlorometanu (B) i temperatury (C) 
w procesie desorpcji, k = 1 J/(m*:s:K), D=1m, T,=413K 
Vai = 3,25:107% kmol/kmol,, Y,,2 = 1,56: 107% kmol/kmol, 
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Ryc. 6. Krzywe wyjścia stężenia n-butanolu (A), tetrachlorometanu (B) i temperatury (C) w procesie 
desorpcji dla różnych temperatur gazu na wlocie do złoża, k=1 J/(m*'s:K, D=1m 


Ya = 3,25:107* kmol/kmol,, Y,.2 = 1,56: 10 * kmol/kmol, 


Na krzywych wyjścia n-butanolu uzyskanych dla temperatury 393 i 413K 
widoczne są dwa maksima stężenia. 

Wpływ średnicy kolumny na przebieg procesu desorpcji, przy stałej 
wartości współczynnika przenikania ciepła k, przedstawiono na rycinie 7. 
W obszarze małych średnic kolumny obserwuje się duże zmiany temperatury 
gazu na wylocie z kolumny ze zmianą jej średnicy. Przy dużych średnicach 
złoża wpływ ten jest niewielki. Średnica kolumny adsorpcyjnej wywiera 
również wpływ na wartość maksymalnych stężeń rozpuszczalników na wylocie 
z kolumny adsorpcyjnej. 
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Ryc. 7. Krzywe wyjścia stężenia n-butanolu (A), tetrachlorometanu (B) i temperatury (C) w procesie 
desorpcji dla różnych średnic kolumny adsorpcyjnej, k = 1 J/(m*'s:K), 7, = 393 K 


a= 325105" kmol/kmol,, 1442 = 1,56: 107% kmol/kmol, 


W laboratoryjnych kolumnach adsorpcyjnych o małej średnicy intensyw- 
ność wymiany ciepła przez ściankę kolumny odgrywa bardzo ważną rolę [2]. 
W celu zbadania wpływu strat ciepła na przebieg desorpcji w kolumnie 
przemysłowej wykonano obliczenia dla dwóch różnych wartości współczyn- 
nika przenikania ciepła przy średnicy kolumny adsorpcyjnej D= 10 m. 
Przyjęto następujące wartości współczynnika przenikania ciepła: k = O (kolum- 
na adiabatyczna) oraz k = 10 J/(m*:s:K). Wyniki obliczeń przedstawiono na 
rycinie 8. Jak widać, przy dużych średnicach kolumny intensywność przepływu 
ciepła przez Ściankę kolumny jedynie w niewielkim stopniu wpływa na 


8 — Szczecińskie Roczniki Naukowe t. XI, z. 1 
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Ryc. 8. Krzywe wyjścia stężenia n-butanolu (A), tetrachlorometanu (B) i temperatury (C) w procesie 
desorpcji dla różnych wartości współczynnika przenikania ciepła k, D=1m, 1, = 393K, 


Yar = 3,25:107* kmol/kmol,, Y;. = 1,56: 107 * kmol/kmol, 


maksymalne wartości stężeń na wylocie z kolumny. W przypadku kolumn 
przemysłowych o dużej średnicy uzasadnione jest więc przyjęcie założenia, że 
adsorber pracuje w warunkach adiabatycznych. 


4. WNIOSKI 


Analiza teoretyczna procesu adsorpcji mieszaniny dwuskładnikowej n-bu- 
tanol-tetrachlorometan wykazała, że lepiej adsorbującym się związkiem jest 
n-butanol, który wypiera tetrachlorometan ze złoża węgla aktywnego. 

Krzywe wyjścia stężenia w procesie desorpcji mieszaniny rozpuszczalników 


LIŻ 


różnią się od krzywych wyjścia dla czystych rozpuszczalników; brak jest 
wyraźnie rozwiniętego obszaru plateau, który obserwuje się w czasie desorpcji 
pojedynczego składnika. 

Temperatura gazu na wlocie do złoża oraz Średnica kolumny adsorpcyjnej 
wywierają wpływ na przebieg desorpcji. Ze wzrostem temperatury gazu na 
wlocie do złoża wzrasta wartość maksymalnych stężeń rozpuszczalników na 
wylocie z kolumny oraz zmniejsza się czas desorpcji. Przy dużych średnicach 
kolumny adsorpcyjnej wartość współczynnika przenikania ciepła od złoża do 
otoczenia jedynie w niewielkim stopniu wpływa na przebieg procesu desorpcji. 
Przy małych średnicach kolumny zmiana średnicy wywiera znaczny wpływ na 
przebieg procesu desorpcji. Przy dużych średnicach wpływ ten jest niewielki. 


5. SPIS OZNACZEŃ 


a,  — stężenie składnika j w fazie stałej [kmol/kg] ([kg/kg], równanie 17) 
stężenie składnika j w fazie stałej odpowiadające całkowitemu po- 
kryciu powierzchni jednocząsteczkową warstwą składnika j [kg/kg] 


R | 


C.;, - ciepło właściwe składnika j w fazie ciekłej [kJ/(mol:K)] 
C,, ciepło właściwe składnika j w fazie gazowej [kJ/(kmol:K)] 
2 
C,=Cat 2 1;,C,, — ciepło właściwe fazy gazowej [kJ/(kmol,:K)] 
j=1 
C,  — ciepło właściwe gazu inertnego [kJ/(kmd,;:K)] 
2 
C,=C,+ Q a,C,, — ciepło właściwe fazy stałej [kJ/(kg: K)] 
g=H 
C,  — ciepło właściwe adsorbentu [kJ/(kg:K)] 
D  — średnica złoża [m] 
G  — gęstość strumienia molowego gazu inertnego [kmol,/(m* :s)] 
AH, — izosteryczne ciepło adsorpcji (desorpcji) składnika j [kJ/kmol] 
k — współczynnik przenikania ciepła od złoża do otoczenia [kJ/(m*: 
's- K)] 
k,  — stała w równaniu 17 
L  — wysokość złoża adsorbentu [m] 
p;  — prężność par składnika j [Pa] 
t — czas [S] 
T  — temperatura [K] 
1, — temperatura fazy gazowej na wlocie do kolumny w procesie adsorpcji 
[K] 
Ilja — temperatura fazy gazowej na wlocie do kolumny w procesie desorpcji 
[K] 
1, — temperatura otoczenia [K] 
R — współrzędna osiowa [m] 
Y, — stosunek molowy składnika j w fazie gazowej [kmol/kmol,] 
Yaj — Wartość Y, na wlocie do kolumny w procesie adsorpcji [kmol/kmol, | 
Y,aj — Wartość Y;, na wlocie do kolumny w procesie desorpcji [kmol/kmol, | 
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z=x/L — wysokość bezwymiarowa 


00 — gęstość gazu inertnego na wlocie do kolumny [kmol,/m*] 
0, — gęstość nasypowa złoża [kg/m'] 
T = Gt/(00L) — czas bezwymiarowy 


BOGDAN AMBROŻEK 


MODELLING OF ADSORPTION AND DESORPTION PROCESSES 
OF BINARY MIXTURE FROM A MOTIONLESS ADSORBENT 
DEPOSIT 


Summary 


In the paper the results of theoretical analysis of adsorption and desorption processes of 
n-butanol and tetrachloromethane vapour mixture in a motionless active carbon deposit were 
presented. It was assumed that adsorption was carried out from a stream of air, whereas the active 
carbon deposit regeneration was made with a stream of heated nitrogen in an open cycle. The 
research was carried out on the basis of the equilibrium sorption theory which assumes the 
existence of local adsorption equilibrium in every point of a deposit. The adsorption equilibrium 
was descreibed by an extended Langmuir's equation. 

The model of the analysed process constitutes a system of three partial differential equations, 
equations of weight balance for both ingredients of the mixture, and energy balance equations. 
Numerical solutions were obtained by means of the line method. 

The influence of adsorption and desorption processes parameters on the course of time 
changes of adsorbate concentration and temperature at the exit from the column was analysed. The 
influence of solvent concentration in the air at the inlet to the column on the course of the 
adsorption process was discussed, as well as the influence of nitrogen temperature at the inlet to the 
deposit, the diameter of the adsorption column, and the intensity of heat losses throug the wall of 
the column on the output curves of desorption. 

The theoretical analysis of adsorption process of binary mixture of n-butanol and tetra- 
chloromethane showed that n-butanol is a better adsorbing compound which displaces tetra- 
chloromethane from an active carbon deposit. 

Along which the increase in the temperature of the gaseous phase at the deposit entry in the 
process of desorption the maximum value of solvent concentration at the outlet from the adsorber 
increases and the time of desorption decreases. For the high values of an adsorption column 
diameter, which is characteristic for industrial columns, the intensity of heat flow through the 
column wall influences the value of maximum concentration at the outlet of the column only in 
a small extent. 


5BOTĄAH AMBPOXEK 


MOJIEJIMPOBAHNE IIPOLIECCOB AJHCOPBIIMM M |ECOPELMMA 
NBYKOMIIOHEHTHOH CMECH M3 HEIIOJ]BAKHOTO 
CJIOA AJICOPBEHTA 


Pe3oMe 


B paóoTe npercTaBJIeHbi pE3YJIbTATbI TEOPETHUECKOTO AaHaJIH3A IpPOLIECCOB AACOPÓRNKU 
MH AecopÓNUM CMeCH IApoB n-OyTaHOJla A TETpaxJIOpMETaHa B HEIOJHBH'KHOM CJIOe AKTHUBHOTO 
yrjiepona. IIpeaAnoJarajlocb, uro aącopÓNMA NpoHCXOJKIIA H3 IIOTOKA BOZJIYXA, a pereHepalna 


7 


CJIOA AKTHUBHOTO YyTJIepONA IIOTOKOM HaTpeTOTO a3OTA B OTKPBITOM HHKJie. ICCJIE1OBAHHA 
UPOBOJKJIACb Ha OCHOBe TEOpHH YpABHOBEIIEHHOŃ COPÓLNAK, IpeANOJararolień CyliecTBOBaHHe 
MECTHOTO a1COPÓNMOHHOTO paBHOBECHA B KAXKIOŃ TOUKE CJIOA. PaBHOBECHE AHCODPÓLMM OIIACAHO 
C IIOMOILIIbIO paCIIIHDeHHOTO YpABHeHHA JleHTMIOpa. 

MONEJIb AaHaJIHBHUPOBAHHOTO IIpOLECCA ABJIAETCA CHCTEMOŃ Tpex nuhdhepeH1NuaJlbHbIX ypaB- 
HeHHf C UaCTHBIMH IIDOH3BOJIHBIMH; ypaBHeHHi MACCOBOTO OaJIaHca AJIA HBYX KOMIIOHEHTOB CMECH 
A ypaBHeHH4A JHEpreTHUECKOTO ÓaJlaHca. UHCJIEHHO€ pelIIeHHE MOJIEJIM IIOJIY4EHO METOJIOM JIHHHH. 

IIpoaHajrmzupoBaHO BJIMAHHe IapaMeTpoB IDpOLECCOB AHCOPÓNHM H NECOPÓNMM B XOJE 
BpeMeHHbIX H3MEHCHHŃ KOHHEHTPALIMH AHCOPÓATOB MH TeMIIepaTypbl IDM BBIXOĄE H3 KOJIOHHBI. 
PaccMOoTpeHo BJIMAHMe KOHIHEHTPALHH pACTBODPHTEJIEŃ B BO3JIYXE IIDH BXOJĄ|E B KOJIOHHy Ha XOJ 
Hpoliecca a1COpóÓNuM, a Tak»Ke BJIAAHHe TeMIlepaTypbl Aa30TA IpH BXOJE B CJIOH, NUaMeTpa 
a1copÓNuOHHOŃ KOJIOHHbI H AHTEHCHBHOCTH IIOTEpH TETIJIA Uepe3 CTEHKH KOJIOHHBI Ha HCXOJIHbie 
KDpUBblie NECOPÓLIMA. 

TeopeTuueckuń aHaJIu3 Hpolecca aqCOPÓLNKM NBYKOMIIOHEHTHOŃ CMEcH n-OyTaHOJI-TeTpa- 
XJIOpMETAH IIOKA3AJI, 4TO ÓoJlee a1COpÓHPYIOLNHM ABJIAETCA COEJHAMHCHHE h-OyTAHOJI, KOTODPBIŃ 
BbITECHA€T TEeTpaXJIODMETAH M3 AKTHBHOTO YyTJIepoJla. 

C nIoBbIIiieHHeM TeMIlepaTypbi Ta30B0oŃń (ba3bI IpH BXONE B CJIOŃ B IIpolecce recopÓHUM 
yBEJIMUMBAETCA 3HAaUC€HHE MAKCHMAJIBHBIX KOHIIEHTpAHMHi paCTBOPKTEJIEŃ IHIpH BbIXOJE H3 AHCOp- 
6epa u cokpaniaeTca BpeMa recopÓnuui. JIJIa OOJIbIIIHX NAaMETPOB ANCOPÓLNMOHHOŃ KOJIOHHBI, 
4TO XapakTepHO JIJIA IIDOMBIIIIJIEHHOCTHM, HKHTEHCHBHOCTb TEIJIA, HPOXOJNAIIETO Uepe3 CTEHKH 
KOJIOHHBI, TOJIbKO B HEKOTODOŃ CTEIIEHH BJIMAET Ha 3HaU€HHE MAKCHMAJIbHBIX KOHIIEHTPALIMŃ IHIPpH 
BbIXOJĄE H3 KOJIOHHBI. 
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MODELE ROZPRZESTRZENIANIA ZANIECZYSZCZEŃ 
W ATMOSFERZE 


Instytut Inżynierii Chemicznej i Chemii Fizycznej Politechniki Szczecińskiej, 
al. Piastów 42, 71-065 Szczecin 


Streszczenie: Substancje toksyczne do atmosfery są emitowane przez źródła 
punktowe, liniowe bądź powierzchniowe. Rozprzestrzenianie zanieczyszczeń w skali 
regionalnej lub mniejszej zwykle opisuje się stosując modele dyfuzyjne. Stanowią je 
cząstkowe równania różniczkowe typu parabolicznego albo eliptycznego. W pracy, 
na podstawie teorii tych równań [6], opisano propagację zanieczyszczeń ze źródeł 
punktowych. Źródła takie mogą działać chwilowo albo ciągle. Rozpatrzono oba 
przypadki. 


Słowa kluczowe: ochrona atmosfery — rozprzestrzenianie się zanieczysz- 
czeń — modele dyfuzyjne — związki fluoru. 


1. ŹRÓDŁA CHWILOWE 


Rozprzestrzenianie zanieczyszczeń emitowanych ze źródeł chwilowych 
zwykle opisuje się równaniami typu parabolicznego o postaci: 


OS. 0 OS| 0 OS| 0 ZŃ 


UE dy 
gdzie 
S — stężenie substancji w powietrzu [g/m*], 
t — czas [S], 
x, y,z — współrzędne układu kartezjańskiego [m], 


f(x, y,z,t) — gęstość źródeł masy [g/m**s], 
D(x), D(y), D(z) — współczynniki dyfuzji [m*/s]. 

Jeżeli ośrodek jest jednorodny, to D(x) = D(y) = D(z) = D. W przypadku 
różnych stałych współczynników dla współrzędnych oznaczamy je przez 
D,, D,, Da. Jeśli gęstość źródeł masy jest równa zeru, otrzymujemy równanie 
jednorodne: 


S, = 04 8,, +02 S,, +03 S,,, (2) 


w. którym qq =D,, 4; = D,.sdg:=D.. 
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Aby uzyskać jednoznaczne rozwiązanie równania dyfuzji, należy określić 
warunki graniczne. Dla równania parabolicznego warunek początkowy polega 
tylko na podaniu wartości funkcji S(x, y,z,t) w chwili t=0. Warunki 
brzegowe mogą być różne, zależnie od rozkładu stężenia na brzegach obszaru. 
Najczęściej są nimi: 

— stężenie na granicy obszaru L-(t), 


m : OS 
— wartość pochodnej, np. 37” 
z 


— związek liniowy pomiędzy pochodną i funkcją. 

Jeżeli znamy równanie różniczkowe i warunki graniczne, możemy znaleźć 
rozwiązanie S(x, y, z, t). Dla skrócenia zapisu rozważymy zagadnienie z jedną 
zmienną przestrzenną. W tym przypadku 


S, = aj S,., +/(x, t). (3) 
Przyjmujemy warunek początkowy w formie 
S(x,0)=0, (4) 
a warunki brzegowe 
S (0, t) =u(l, t) =0. (5) 


Rozwiązanie tego zagadnienia dla odcinka I i czasu działania źródła t ma 
postać: 


S(x, 1) = IG y ap - (7) 2-9 |sn Tesi z c/e. 2) dt dt = 
n=l 
00 


tl 


= | |se. ć, t—a) fl, o)dłdz, (6) 


00 


gdzie © oznacza położenie źródła emisji na współrzędnej x. 

Wyrażenie G(x, t, 6) nazywamy funkcją źródła. Funkcja źródła, rozpat- 
rywana jako funkcja x, przedstawia rozkład stężenia dyfundującego składnika 
na odcinku O < x < l w chwili t, jeżeli stężenie w chwili początkowej równa się 
zeru i w tej chwili w punkcie x=© wydziela się jednostkowa masa tego 
składnika, a na końcach odcinka podtrzymywane jest przez cały czas stężenie 
zerowe. Należy zauważyć, że dla wszystkich O<x<li t>0 funkcja 
G(x, t, 6) >0. 

Jeżeli moc źródła wynosi Es, to 


S(x, t) = Eg G(x, t,), (7) 


co oznacza, że rozwiązanie zagadnienia dyfuzji polega na znalezieniu funkcji 
źródła G(x, t, £). 
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W przypadku rozprzestrzeniania się zanieczyszczeń w atmosferze dla 
współrzędnej x interesuje nas funkcja źródła dla prostej nieskończonej. Jest ona 
określona wzorem: 


AMIE! _(r—03 
G(x, t, 6) = r 7» zza | (8) 


Zależność (8), wyrażoną funkcją Gaussa, można nazwać rozwiązaniem pod- 
stawowym równania dyfuzji. 
Funkcję źródła dla współrzędnej y, przez analogię, zapisujemy w postaci: 


l (yn) 
G 3 2 = —— —— + 2 
(y t 9) , +/możtt €Xp | 2 (9) 


gdzie yn jest położeniem źródła emisji na współrzędnej y. 

Nieco inaczej konstruuje się funkcję źródła dla współrzędnej z (wysokości). 
Współrzędna z zmienia się w zakresie 0 < z < w. Zagadnienie sprowadza się 
zatem do poszukiwania rozwiązania równania dyfuzji dla półprostej. Ze 
względu na to, że efektywne wysokości źródeł emisji są relatywnie małe 
w stosunku do zakresu zmian współrzędnej z, trzeba rozpatrywać przypadek, 
gdy źródło leży blisko początku półprostej. Prowadzi to do poszukiwania 
rozwiązania równania dyfuzji dla t > 0, przy warunku początkowym 


S(z, 0) = p(z) (10) 


i zależnie od charakteru procesu jednego z trzech warunków brzegowych. 
W teorii propagacji zanieczyszczeń w atmosferze zwykle rozpatruje się drugi 
warunek brzegowy w postaci 0S(0, t)/0z = O. 

Funkcję źródła przy jednorodnym warunku brzegowym drugiego rodzaju 
i warunku początkowym (10) opisuje się wzorem: 


z l (z-H)? (z+H)? 
G(z, Bzz (e | - dażt |+er| - dażt |; (11) 


Przez H oznaczono wysokość efektywną emitora. Sposób jej obliczania został 
podany w pracy [5]. 
Z teorii równań parabolicznych [6] wynika, że funkcję źródła dla półprze- 
strzeni można określić według zależności: 
LPĄ2 
A Ę śe 


2 
dar t 


1 3 
G(6,92,64,1, HB) = 6,6, 6, = (zz) 
ZR BRE wj 


OLq OL2 Ol3 


zexp| Uh ap| -5-7 +ap| Gp (12) 


4aż t 4a3 t 4ażt 
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Równanie (12) jest słuszne, jeżeli położenie środka wyemitowanej chmury 
w czasie rozprzestrzeniania jest stałe (ośrodek nieruchomy). Jeśli masy po- 
wietrza przemieszczają się w czasie rozprzestrzeniania, np. wzdłuż współrzędnej 
x ze średnią prędkością |9,|, to położenie środka chmury zmienia się według 
zależności 


0=|9,]t. (13) 


Współczynniki «,, «,, «a należy wyznaczyć doświadczalnie. 


2. DYFUZJA W STRUMIENIU I ŹRÓDŁA CIĄGŁE 


Rozpatrzmy dyfuzję jednego składnika w strumieniu niestacjonarnym bez 
źródeł. Dla tego przypadku mamy następujące równanie dyfuzji: 


OS 0 08 ) OŃ) 5) GR) 
a-z|POż|+3|?P0%|+520 |- 

0 0 0 

a 090 0, 8297030). (14) 


Jeżeli przyjmiemy, że w półprzestrzeni z > O dany jest strumień skierowany 
wzdłuż współrzędnej x, który porusza się ze średnią prędkością |9,|, to równa- 
nie (14) upraszcza się do postaci: 


OŚ. 0 0S O ÓS| O OS OŚ 
= |? w |+z|D W |z| oz |-Bl. (15) 


Rozwiązania równania (15) poszukujemy konstruując funkcję źródła 
G = G,G,G,. Funkcję źródła G, buduje się rozpatrując zagadnienie dla prostej 
x, zakładając D(x) = D,. Równanie dyfuzji dla tego przypadku zapisujemy 
następująco: 


08 „628 ,, ,68 


„5 => 04 =D,-0<x< w. (16) 


Podstawiając 0 =ait i B=|9,|/aq uzyskujemy 


08 _0*5_ 08 5 
00 0x © 0x 
Następnie podstawiamy 
u = pe'">, (18) 


obliczamy pochodne 05/00, 0S/0x i 07S/0x*, wynik wstawiamy do (18) 
i otrzymujemy 
óv 0%v 0? 8 


1 WE 1 W, 2 = w. 
3 aż ąz © G 
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gdzie 


19, 
=P 
SB - (zz) 


Równanie (19) rozwiązujemy podstawiając 
v=fe", (20) 


obliczamy pochodne 0v/00 i 0*v/0x*, wprowadzamy je do równania (19) 
i otrzymujemy 


06 06 
60 0x7 SP 
albo 
06 076 
Funkcja źródła dla równania (21) ma postać: 
za: xp] " (22) 
"  24/nażt 4aq t | 


Funkcję źródła dla równania (17) zapisujemy w formie: 


> 1 3” 19, | 
G, = 6,0 -—zap| -gr-z( 2) |. 23 
2 mat | dat dla, ER) 


a funkcję źródła dla równania (16) 


1 (x —|9,| 4) 
G,G=G,xef" =——<=0%p |- (24) 
2.,/najt 4a t 
Funkcję źródła G, konstruuje się na podstawie równania: 
OS OS? 
Ma ona postać 
I A 
G, =— exp | ————— |, (26) 
are, / no, 46) 


2 


[9,| 

Funkcja źródła G, określa się na podstawie podobnego do (25) równania, 
zapisanego dla współrzędnej z. Rozwiązuje się je dla z>0, przy drugim 
warunku brzegowym 0(0, t)/0z = 0. W wyniku uzyskuje się 


=-_1 Śgiak ESR 
aezżpeawiolok © 


D; 
JĆ jest współczynnikiem. 


gdzie 6; = — x jest współczynnikiem. 


=2 
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Współczynniki 6,, 6, są wyrażane w metrach. Wyznacza się je na podstawie 
danych doświadczalnych. 
Funkcję źródła dla przestrzeni trójwymiarowej uzyskuje się z zależności 


G(x, ), z, t; 8,, 9, H) = G,G,G,. (28) 


Znając funkcję źródła możemy zapisać równanie określające czasowo-prze- 
strzenny rozkład stężenia w niestacjonarnym strumieniu powietrza: 


Jes = gE (x —|[9,| 2)? | (7-1) 
Ś > 43 t)=|——= == 2 Pu” = r a 
GZVAZY) Z) zazjo| 462 PB 40; 


z-H)? z+H)? 
em - z |+ep|- Gz | (29) 
gdzie 
6, =ojt [m], 
E — ilość wyemitowanej substancji [kg]. 


Równanie opisujące rozkład przestrzenny stężenia w strumieniu dla punk- 
towego źródła ciągłego (t > co) o stałej wydajności g(t) = E otrzymuje się 
z rozwiązania niejednorodnego równania dyfuzji następująco: 


[e 0) 


S(x, y, z) = je. G,G,o(t)dt = EG, G, |G.d — 


0 


0 
= AS 240 | gi (z-H)? (+H)? 
aero]. © 


gdzie E oznacza strumień masy emitowanej substancji [kg/s]. W obliczeniach 
zwykle przyjmuje się je jako stężenie okresowe, chwilowe i średnie roczne. 
Wtedy E oznacza średnie natężenie emisji w ciągu np. 30 min lub roku. 

Zależności (29) i (30) uzyskano przy założeniu, że emitowane zanieczysz- 
czenia nie ulegają przemianom w atmosferze, masy powietrza poruszają się bez 
zakłóceń wzdłuż współrzędnej x ze średnią prędkością |9,|, współrzędne 
zmieniają się w zakresie: —0o <x<o«, —w<y<o, z=20 i dla z jest 
spełniony drugi warunek brzegowy. Są więc one słuszne tylko dla tych 
warunków. 

Przedstawiony opis, pod względem formalnym, można zaliczyć do modeli 
Eulera [9]. W szczególności jest propozycją prostego sposobu uzyskiwania 
rozwiązań analitycznych. Wielkości 6,, 6, spełniają rolę współczynników 
modelu. Zwykle opisuje się je równaniami empirycznymi [7, 9]. 

Interesujące jest porównanie otrzymanych zależności z powszechnie sto- 
sowanym modelem Gaussa. Z matematycznego zapisu końcowych równań 
obu modeli wynika, że stosując je uzyska się te same rezultaty podczas obli- 
czeń przy spełnieniu pomiędzy wartościami liczbowymi tych współczynników 
zależności: 


12) 


0, = „/26,, (31) 
0, =,/26,. (32) 


Podczas weryfikacji wystarcza zatem sprawdzić jeden z tych modeli, ponieważ 
drugi, po uwzględnieniu zależności (31) i (32), będzie także spełniany. 


3. SYSTEMY OCENY JAKOŚCI POWIETRZA 


Według obowiązujących norm jakość powietrza ocenia się na podstawie stę- 
żeń chwilowych, średnich dobowych i średnich rocznych [2, 4]. Zwykle określa 
się tylko stężenia chwilowe i średnie roczne [1]. Stężenia te można wyznaczać 
doświadczalnie albo obliczyć stosując odpowiednie modele matematyczne. 

Obliczając rozkłady stężeń wykorzystuje się odpowiednie systemy. Zawie- 
rają one równania opisujące czasowo-przestrzenne rozkłady stężeń, zależności 
korelacyjne określające parametry tych równań oraz odpowiedni sposób 
uwzględniania warunków meteorologicznych. Znane są liczne systemy, np. 
Julicha, Pasquille'a, Brigsa, St. Luis, Kluga [8]. Parametry o,, 0,, o, oblicza się 
z równań empirycznych. Parametry o,, o, zwykle opisuje się wzorem: 


0,=0,=Ax", (33) 
gdzie 
X — odległość od źródła [m], 
A,a — współczynniki empiryczne. 


Więcej typów zależności stosuje się do opisu parametru o, [8]. Do bardziej 
znanych należą równania: 
Pasquille'a, Gifforda, Singlera, Smitha (PGSS) oraz Kluga 


G=dax, (34) 
Martina-Tikvarta 
0, =ax+C, (35) 
Brigsa 
0, =ax(1+bx), (36) 
gdzie 
x — odległość od źródła [m], 
a,b,c — współczynniki empiryczne. 


Obowiązujące w Polsce wytyczne [1] oparte zostały na systemie Pasquille'a 
wykorzystującym model Gaussa. Stosowane w tym systemie równanie, opisu- 
jące rozkłady stężeń, ma postać: 


[e 0) 


E (Wn) 
S(x, y, z) = | S(x,y, z, tydr=—— mh 
ŻY, | GO zawie? | 2aż | 


« |ew| --7 +xp| z [g/m]. (37) 


0 


5) ż , 
207 207 


126 


Parametry równania (37) w wytycznych opisano następująco: 
0,=0,=Ax*, (38) 
0,=Bx*, (39) 


gdzie 

a = 0,367 (2,5 — m), 

b = 1,55exp (— 2,35 m), 

A =0,08[6m **+1—1ln(H/z5)], 

B = 0,38 m'** [8,7 — ln (H/z5)], 

m — współczynnik meteorologiczny (dane tabelaryczne), 

H — wysokość efektywna emitora, 

z — współczynnik korekcyjny dla terenu (dane tabelaryczne). 

Uważamy, że przed przystąpieniem do obliczeń stężeń imisyjnych trzeba 
dla danego regionu zweryfikować metodykę przez porównanie stężeń ob- 
liczeniowych i zmierzonych. Za substancję modelową należy przyjmować taką, 
która jest emitowana tylko z jednego obiektu; nie powinna występować w tle 
i nie może być przenoszona z innych regionów, ponieważ udziały te są trudne 
do oszacowania. Trzeba zauważyć, że zgodność wartości stężeń obliczonych 
z modelu matematycznego i stężeń zmierzonych zależy nie tylko od jakości 
modelu matematycznego, lecz także precyzji pomiarów wykonywanych na 
poszczególnych punktach sieci monitoringu oraz dokładności w określeniu 
wielkości emisji. 

Takie obliczenia sprawdzające przeprowadzono dla rejonu Szczecina w la- 
tach 1992— 1993. Za substancję modelową przyjęto rozpuszczalne w wodzie 
związki fluoru, ponieważ w regionie są one emitowane głównie przez jeden 
zakład. Rozkłady stężeń chwilowych S$ i średnich rocznych S$, związków 
fluoru w atmosferze obliczono na podstawie wyznaczonych doświadczalnie 
wielkości emisji i parametrów technicznych emitorów. Dla przykładu w tabeli 1 
podano takie dane dla 1993 r. 


Tabela 1. Charakterystyka źródeł emisji związków fluoru 


owa ż LIES Ś Ę z a: 
emitora [m] [kg/h] [m] [m/s] [K] [m] | 
Po-1 0 0,0645 1,20 10,5 333 45 
Po-2 80 0,0657 0,60 195 313 25 
Po-3 145 0,1200 1,20 10,5 333 45 
Po-4 165 0,1363 0,60 19,5 313 25 
Po-5 10 0,8300 0,31 5,8 311 i 
Po-6 265 0,3404 1,60 13,8 340 45 
Po-7 290 0,0680 0,64 12,7 325 32 
Po-8 —20 0,2670 3,96 112 328 60 
Po-9 24.5 2,7412 3,96 12,4 328 60 
Po-10 —575 0,1524 3,20 10,2 303 42 
Po-11 —575 0,2836 4,20 6,0 318 42 
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Tabela 2. Porównanie stężeń obliczonych i zmierzonych rozpuszczalnych w wodzie związków 
fluoru. Rejon Polic 1992 r. 


Miejscowość p FE" | osSp/% 
S;o[Hg/m ] | S.[Ug/m'] | S„[Hg/m'] 


Jasienica 
Tatynia 


Police — Fabryczna 
Police — Stadion 
Police — Hotel 
Wieńkowo 


S30 — obliczone stężenie chwilowe, S$, — obliczone stężenie średnie roczne, S,, — zmierzone 
stężenie średnie roczne. 


Tabela 3. Porównanie stężeń obliczonych i zmierzonych rozpuszczalnych w wodzie związków 


fluoru. Rejon Polic 1993 r. 
SzplSa | 
1,4 


EEE Pomiar 
Miejscowość ZUJ PRRZĄTYTET 
S30 EE EZ [ug/m ] | S„[Hg/m ] 


Jasienica 0,07 
Tatynia 

Police — Fabryczna 

Police — Siedlecka 


Police — Hotel 
Wieńkowo 


Objaśnienia jak do tab. 2. 


Na rycinach I i 2 przedstawiono w postaci wybranych izolinii rozkłady 
stężeń Średnich rocznych. Na rycinach zaznaczono także położenie stacji 
pomiarowych. Odczytane z takich rycin, dla poszczególnych stacji, obliczone 
wartości stężeń chwilowych i średnich rocznych podano w tabelach 2 i 3. 

Zmierzone wartości stężenia średnich rocznych w poszczególnych stacjach 
wyznaczono na podstawie 24-godzinnych pomiarów ciągłych. Obliczano 
najpierw wartości średnie dla każdego miesiąca, a potem roku. Wartości te 
podano w tabelach 2 i 3. 

Dla części stacji uzyskano dobrą zgodność, natomiast dla kilku stanowisk, 
różnych w latach 1992 1 1993, rozbieżność była znaczna. 

Mała ilość końcowych wyników, warunkowana niewielką liczbą stacji 
sieci monitoringu, nie daje podstaw do opracowania ich metodami statys- 
tycznymi. 

Średnia wartość stosunku stężeń wyznaczonych z pomiarów i obliczonych, 
wynosząca około 1,8, wskazuje na błąd systematyczny. Może on dotyczyć 
modelu albo sposobu określania stężeń średnich rocznych z pomiaru (dokład- 
ność metody analitycznej, liczba wykonanych oznaczeń itp.). Nie ma prostych 
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ję -3200 -2400 -1600 —800 0 800 1600 2400 3200 400% 
| A 


-3200 —2400 —1600 —800 


44000 -3200 -2400 >1600 =800 o BOQ0 1800 2400 3200 $ 
X[m). 


Ryc. 1. Wybrane izolinie stężeń średnich rocznych S,* 100 [ug/m*] rozpuszczalnych w wodzie 
związków fluoru. Rejon Polic 1992 r. 


-3200 -2400 -1600 —800 0 B00 1600 2400 3200 40B 
A 


34000 -3200 -2400 -1600 —800 (0 BOO 1600 2400 3200 aż 
X[m] 


Ryc. 2. Wybrane izolinie stężeń średnich rocznych S,* 100 [ug/m*] rozpuszczalnych w wodzie 
związków fluoru. Rejon Polic 1993 r. 
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metod, aby to rozstrzygnąć. Sposób pośredni może polegać na jednoczesnym 
badaniu na tym samym obszarze kilku substancji modelowych [3]. Według 
takiej oceny omawiany błąd należy przypisać metodzie wyznaczania średnich 
rocznych stężeń związków fluoru z pomiarów. 


4. WNIOSKI 


Spośród stosowanych przy opisie propagacji zanieczyszczeń w atmosferze 
modeli do podstawowych należą: różniczkowe modele Eulera, całkowe modele 
Lagrange'a oraz modele Gaussa. 

Stosując modele różniczkowe w postaci cząstkowych równań różniczkowych 
parabolicznych ieliptycznych opisano rozprzestrzenianie się zanieczyszczeń emi- 
towanych z punktowych źródeł chwilowych i ciągłych, a w szczególności przed- 
stawiono prosty sposób uzyskiwania, ważnych dla rozważań teoretycznych, 
rozwiązań analitycznych. Polega on na zastosowaniu pojęcia funkcji źródła. 

System obliczeń sprawdzano w latach 1992 —1993, porównując obliczone 
i wyznaczone z pomiarów na stanowiskach istniejącej w rejonie Polic sieci 
monitoringu stężenia Średnie roczne rozpuszczalnych w wodzie związków 
fluoru. 


JERZY STRASZKO, LUCIANO A. PAULO, ARTUR WITKOWSKI 
MODELS OF POLLUTION PROPAGATION IN THE ATMOSPHERE 


Summary 


The description of the pollution propagation process in the atmospheric air from occasional 
and continuous sources was presented. Propagation models in an analytical form were obtained 
basing on the theory of partial differential equations of the parabolic type, and the utilisation of the 
source function notion in particular. 

The calculation system, in the form of computer program, was obtained by adding the 
empirical formulae characterising transport coefficients to the models and taking meteorological 
conditions into consideration in an appropriate way. The obtained system was tested by 
comparing the average yearly concentration of water-soluble fluorine compounds calculated and 
measured by the network of monitoring system existing around the "Police" Chemical Plant. In the 
calculations all existing emission sources of these compounds in the Szczecin urban area were 
included, i.e. phosphatic fertilisers plants and coal burning processes in power plants and 
thermal-electric power stations. 


EXKM CTPAIIIKO, JTYHMAHO A. IAYJIO, APTYP BATKOBCKH 
MOJIEJIM PACIIPOCTPAHEHHA 3ATPA3HEHHA B ATMOCOEPE 


Pe3oMe 


IIpe4cTaBJieHo ormcaHne HpoHecca pacnHpocTpaHeHia 3aTpA3HeHHi B aTMOCÓEPHOM BO3NYXE, 
U3JIYUAeMPIX H3 HEIIOCTOAHHbIX H IIOCTOAHHbIX ACTOUHHKOB. MOJHEJIM pacHpocTpaHeHHA B AHaJIH- 
TUUECKOM BHJIE IOJIYHCHbl Ha OCHOBe Teopuu NudhdhepH1uaJlbHbiIX ypaBHeHHA C UACTHBIMH 


9 — Szczecińskie Roczniki Naukowe t XI, z. 1 
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HpOH3BOJHBIMA IlapaOOJIMUECKOTO THIIA, A B UACTHOM CJIYUA€ OBIJIO IIDUMEHEHO IIOHATye (DyHKIIAA 
HACTOUHKKA. 

CucTeMbl BbIUKCJIEHHŃ, B BHUJIE KOMIIBIOTEPHOŃ IPOTPAMMPBI, IIOJIyUEHbI IIyYTEM IIDHCOEJTHHE- 
HHA K MOJIEJIAM 3MIIADHUECKHX (DODMYJI, ONDEJIEJIAIOLIAX KOJÓÓUNAEHTbBI TPAHCIIODTA H YUKTbI- 
BAA COOTBETCTBYFOLIHM OÓpa30M METeOYyCJIOBHA. IIojryueHHaa CHCTEMA ObiJIa IpoBepeHa IIyTeM 
CpaBHeHHA BbIUMCJIEHHbIX M H3MEPpCHHbIX pE3YJIbBTATOB HPA IIOMOINM CylieCTBYFoLiech B palioHe 
XuMHUECKOTO KOMÓWHATa „„Ilojrmiie' ceTA MOHHTODHHTA CpEHHETOHOBBIX KOHIIEHTDANHA pacT- 
BOpHUMBIX B BOJIe COEĄUHeHHMA (bTOpa. IIpH BbIUMCJIEHAAX ObIJIM YUTEHbI BCE, HAXOJALNIAECA 
B IllerruHckoń arJJoMepanuu, HCTOUHHKH H3JIYUCHHA 3THUX COENAHCHHŃ, TAK HA3bIBAeMblie (hPaÓpHKH 
(bochopHBbix ynoOpeHuń, a TaK?KE IIDOIIECCI CYKHATAHHA YTJIA B JJIEKTDOCTAHIIMAX H TEIIIOJJIEKTDOC- 
TAHNHAX. 
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